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1.0  Introduction 

The  Development  and  Evaluation  of  Fusion  Techniques  (DEFT)  in-house  program  was  initiated  with 
the  intent  of  addressing  multiple  areas.  The  first  purpose  of  DEFT  was  to  extend  the  development  of 
algorithmic  concepts  that  had  been  identified  and  explored  under  the  previous  Complementary  Advanced 
Fusion  Exploration  (CAFE)  in-house  program.  CAFE  investigated  promising  recent  innovations, 
primarily  involving  mathematical  concepts  within  the  area  of  basic  research,  and  did  a  preliminary 
assessment  of  their  validity  as  well  as  a  preliminary  qualitative  assessment  of  each  concept’s  potential  for 
application  to  Air  Force  Command,  Control,  Intelligence,  Reconnaissance,  and  Surveillance  (C2ISR) 
capabilities  and  systems.  DEFT  was  intended  to  continue  the  investigation  of  potential  exploratory 
development  applications  of  the  CAFE  algorithmic  concepts  and  to  investigate  concepts  that  were 
revealed  as  extensions  of  the  original  CAFE  concepts. 

Areas  of  investigation  under  CAFE  that  were  originally  identified  for  further  exploration  under 
DEFT  included  the  use  of  Trilinear  tensors  (including  Homography)  for  use  in  the  reprojection  of  images, 
the  Nash  (Grocholsky)  Equilibrium  Approach  to  Sensor  Management,  the  Virtanen  Methodology 
Approach,  improvements  in  handling  the  Out  of  Sequence  Problem  (OOSP)  and  Out  of  Sequence 
Measurement  (OOSM)  updating,  and  further  study  in  the  area  of  Particle  Filters  and  their  applications. 

In  addition,  DEFT  was  intended  to  study  other  Exploratory  Development  topics  not  previously 
investigated  under  CAFE.  These  topics  included  Bayesian  Network  exploitation  for  fusion,  Optimization 
of  ISR  platform  route  determination,  and  Radio  Frequency  (RF)  tag  utilization. 

The  first  area  to  be  pursued  under  DEFT  was  the  area  of  Particle  Filtering.  As  work  proceeded 
on  this  area,  it  became  clear  that  study  in  this  area  should  become  the  focus  of  the  DEFT  in-house 
program.  This  is  what  took  place.  As  a  result,  no  work  took  place  on  the  other  originally  intended  topics. 

As  the  work  focused  on  Particle  Filtering,  the  directions  of  the  work  were  clarified.  The  nature  of 
the  problem  being  addressed  was  that  traditional  trackers,  based  around  the  Kalman  Filter,  were  designed 
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to  deal  with  linear  targets  with  Gaussian  noise  statistics.  However,  targets  are  not  constrained  to  behave 
in  these  ways,  and  so  methods  are  needed  to  deal  with  targets  that  are  not  behaving  in  these  ways. 

The  Particle  Filter  is  one  such  method.  It  does  not  make  the  assumptions  of  linearity  and 
Gaussianity  that  the  Kalman  Filter  does,  but  it  requires  a  large  number  of  computations,  depending  on 
how  parameters  are  set. 

To  assess  the  usefulness  of  the  Particle  Filter,  a  systematic  program  was  initiated  to  investigate  it 
to  assess  its  practicality  and  usefulness.  This  program  involved  a  number  of  areas.  First,  further  study 
was  required  to  understand  the  mathematical  basis  for  why  the  Particle  Filter  provides  improved  filtering 
and  is  a  basis  of  improved  tracking  in  situations  of  interest.  Next,  it  was  desired  to  create  and/or  obtain  a 
simulation  as  a  framework  for  systematic  parametric  testing  of  the  performance  of  the  particle  filter 
versus  other  types  of  filters.  However,  the  goals  of  this  effort,  in  addition  to  those  of  CAFE  were  to  take  a 
more  focused  approach  to  fusion  techniques.  That  is,  DEFT  looked  to  focus  on  a  deep  mathematical 
understanding  of  algorithms  researched  under  CAFE,  as  well  as  to  go  beyond  the  work  done  under  CAFE 
by  creating  in-house  implementations  for  the  testing  of  these  algorithms.  The  team  looked  to  create  an  in- 
house  testbed  that  would  serve  as  a  model  for  other  multi-sensor  fusion  studies,  as  well  as  to,  as  the  name 
suggests,  evaluate  the  performance  of  these  methods.  These  tasks  were  created  to  advance  and  assess  new 
algorithmic  concepts  in  order  to  enhance  the  capabilities  available  to  track  maneuvering  and 
unpredictable  targets.  DEFT  searched  for  algorithms  and  techniques  that  would  optimally  fuse  the 
information  available. 

The  scenarios  for  concern  included  those  situations  that  involve  highly  non-linear  and  non-Gaussian 
characteristics.  While  existing  methods,  such  as  the  Kalman  (KF)  and  Extended  Kalman  Filters  (EKF), 
perform  sufficiently  in  linear  and  Gaussian  cases,  these  algorithms  tend  to  lose  tracks  and  perform 
unreliably  when  faced  with  data  that  exhibits  alternative  behavior.  In  order  to  gain  increased  tracking 
accuracy  and  reliability,  the  team  proposed  the  prospects  of  using  a  Particle  Filtering  (PF)  algorithm. 
Based  on  research  done  under  CAFE,  as  well  as  the  fact  that  there  are  no  assumptions  of  linearity  or 
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Gaussianity  made  in  the  PF,  these  filters  showed  potential  for  providing  an  innovative  approach  to  the 
non-linear  and  non-Gaussian  scenarios  under  which  the  Kalman  Filters  have  shown  possible  weaknesses. 

The  approach  of  this  effort  was  to  mathematically  characterize  the  algorithms  under  investigation.  Once  a 
working  definition  and  understanding  of  the  techniques  was  developed,  the  team,  in  collaboration  with 
Syracuse  University,  created  an  implementation  of  these  algorithms  to  evaluate  the  performance  of  the 
techniques  on  the  standards  of  position  error,  computational  expenditures,  and  characterizations  of 
situations  under  which  a  filter  showed  improved  tracking  results. 

2.0  Kalman  Filtering 

Developed  in  the  1960’s,  the  Kalman  Filter  (KF)  was  the  first  of  a  family  of  algorithms  aimed  at 
estimating  the  current  state  of  an  object.  The  KF  is  an  optimal  (unbiased)  estimator  for  objects  with  linear 
and  Gaussian  characteristics,  for  this  algorithm  assumes  that  the  object  has  these  attributes.  The  basic 
theory  is  to  predict  the  current  state  of  an  object,  upon  reception  of  observations,  by  taking  the  weighted 
sum  of  the  estimate  and  the  most  recent  observations.  Hence,  as  the  general  theory  suggests,  the 
algorithm  is  broken  up  into  prediction  and  update  stages. 

2.1  Algorithm 

Prediction 

The  equations  which  calculate  the  predicted  state  vector  and  covariance  matrix  are: 


x  =  F  ■  x  +  OWN 

(1) 

P  =  F  P -Ft  +Q 

(2) 

o  x  and  x  are  the  state  vector 
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o  P  and  P  are  error  covariance 
o  F  is  the  state  transition  matrix 
o  Q  is  process  noise  covariance 
o  OWN  is  a  control  vector 

Note  that  when  a  matrix  or  vector  symbol  is  covered  with  a  hat,  the  parameter  represents  a  measurement- 
updated  estimate ;  when  a  matrix  or  vector  symbol  is  covered  with  a  tilde,  this  represents  a  parameter  that 
is  a  time-updated  prediction.  Therefore  the  initial  state  estimate  will  be  represented  by  x  and  the  initial 

error  covariance  matrix  will  be  represented  by  P  . 

As  a  recursive  algorithm,  the  process  must  start  with  an  initially  defined  state  estimate  and  error 
covariance.  The  elements  of  the  state  vector  are  the  parameters  that  are  to  be  estimated.  While  the  state 
vector  could  include  a  wide  range  of  components,  when  tracking  ground  targets  the  parameters  of  interest 
are  typically  the  position  and  velocity  of  the  object  in  Cartesian  coordinates.  A  depiction  of  the  state 
vector  may  appear  as  below. 


x 

x 


x  = 


y  > 
y 


(3) 


A  dot  above  a  variable  indicates  the  first  derivative  of  motion  in  that  direction  (also  called  the  velocity  in 
that  direction).  Now,  the  error  covariance  matrix  provides  information  about  the  accuracy  of  the  state 
estimate,  or  can  be  viewed  as  a  measure  of  how  the  components  of  the  state  vector  vary  with  respect  to 
each  other.  In  a  more  mathematically  thorough  statement,  according  to  the  definition  of  state  error 
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covariance,  the  error  covariance  matrix  will  be  an  ordered  set  of  the  covariance  between  each  ordered  pair 
of  differences  among  the  actual  state  and  the  estimated  components  of  the  state  vector. 


P  =  £[(x-x)(x-x)r] 

(4) 

P  =  £[(x-x)(x-x)2’] 

(5) 

Similarly,  the  process  noise  covariance  matrix,  Q  ,  takes  into  account  the  error  involved  in  the  model  that 
is  chosen  for  the  system.  This  error  can  be  incorporated  into  the  system,  by  making  assumptions 
concerning  the  models  chosen  for  the  system,  or  as  will  be  observed  in  the  following  section,  through  an 
approximation  of  the  system  parameters  that  does  not  include  all  higher  order  derivatives. 

The  state  transition  matrix/7  is  essential  to  understanding  the  Kalman  Filter.  This  is  the  matrix  that,  when 
multiplied  by  the  state  vector,  will  produce  the  predicted  state  of  the  object  at  the  subsequent  time  step. 
Since,  this  research  is  primarily  concerned  with  tracking  and  motion,  the  matrix  that  will  be  used  most 
often  in  this  paper  is  easily  derived  from  the  common  kinematic  equations  found  in  physics.  Yet,  it  should 
be  noted  that  the  KF  can  be  applied  to  other  applications,  and  therefore  may  have  other  appropriate 
transition  functions.  The  equations  that  will  be  used  are  as  follows: 


+  x(t_1  ■  dt  +  \/ 2xk_l  ■  dt2 ... 

(5) 

II 

*• 

+ 

X: 

(6) 

The  ellipses  indicate  the  remaining  derivatives  of  motion,  with  which  we  are  not  concerned  due  to  our 
linear  assumption.  Taking  the  coefficients  of  the  elements  of  state  vector,  namely  the  coefficients  of 
velocity  and  position,  the  state  transition  matrix  can  be  derived  from  these  equations  as: 
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(7) 


1  dt  0  0 

0  10  0 

0  0  1  dt 

0  0  0  1 


The  ownship  or  control  vector,  which  adjusts  the  prediction  for  a  sensor  mounted  on  a  moving  platform 
(such  as  on  an  aircraft  or  a  car),  is  optional.  Although  some  refer  to  this  as  the  multiplication  of  a 
transition  matrix  B  and  a  state  vector  u  with  components  of  the  ownship's  velocity  and  acceleration 
which  are  assumed  to  be  known,  here  it  will  be  referred  to  as  OWN.  As  described,  the  OWN  is  often 
represented  as: 


B  u 


dt 

1 


dr_ 

2 

dt 


velocity 

acceleration 


(8) 


Update 

The  update,  upon  the  reception  of  a  measurement,  corrects  the  time  updated  estimate  and  error 
covariance.  These  formulas  are  derived  from  the  minimization  of  the  MSE  of  the  state  and  the 
measurement,  making  the  Kalman  Filter  an  optimal  estimator.  The  update  formulas  are  as  follows: 


K  =  P  Hr[H  P  Hr  +R]  1 

(9) 

P  =  (I  -  K  •  H)P 

(10) 

x  =  x  +  K(z  -  H  •  x) 

(11) 

o  K  is  the  Kalman  gain 

o  P  and  P  are  error  covariance 
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o  H  is  the  state  to  measurement  transformation  matrix 
o  z  is  the  measurement  vector 
o  R  is  measurement  noise  covariance 
o  x  and  x  are  the  state  vector 

Within  the  update  stage,  the  measurement-updated  estimate,  which  is  used  as  an  input  for  the  next 
iteration  of  the  filter,  is  calculated  by  taking  the  sum  of  the  time -updated  estimate  and  a  weighted 
difference  known  as  the  innovation.  The  innovation,  also  called  the  residual,  is  the  difference  between  the 
measurement  z  and  the  product  H  •  x .  The  matrix  H  is  necessary  because  it  transforms  a  vector  from  the 
state-space  to  the  measurement-space.  Otherwise  these  vectors  will  likely  not  be  the  same  size  or  be  in  the 
same  scale.  Given  information  from  N  sensors  centralized  measurement  fusion  is  performed  in  the  KF  as: 

Hk=[HlSi,H{Si,...,HlsJ 

R/c,Sy  “  *  0 

Rk=  :  \  : 

0  Rk,SN 

where  Sx  SN  indicates  the  sensor  from  which  information  is  related  and  Rk  is  the  measurement 
covariance  matrix. 

The  weight  that  is  calculated  is  known  as  the  Kalman  gain,  which  takes  into  account  the  variance  of  error 

in  the  estimate,  P ,  and  the  variance  of  noise  in  the  measurement,  R .  Notice,  as  the  estimate  error 
covariance  goes  to  zero,  the  Kalman  gain  goes  to  zero,  indicating  that  the  estimate  is  highly  reliable.  This 

will  result  in  a  P  with  similarly  smaller  values  as  P  and  a  measurement-updated  estimate  that  is  very 
close  to  the  value  of  the  time -updated  estimate.  Alternatively,  if  the  measurement  noise  covariance 

approaches  zero,  the  Kalman  gain  gradually  progresses  towards  the  value  Hr  .  As  a  result,  the  innovation 
will  be  more  highly  weighted  and  the  measurement-updated  estimate  will  be  farther  away  from  the  time- 
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updated  estimate.  Notice  that  still  in  this  case  where  the  Kalman  gain  is  larger,  the  error  covariance  is  still 
being  minimized.  Hence,  the  general  idea  of  the  Kalman  gain  equation  is  that  it  designates  how  much  the 
innovation  will  change  the  predicted  state. 


2.2  Derivation 

Although,  we  have  preliminary  described  a  few  of  the  parameters  that  will  be  used,  for  this  derivation 
some  of  the  terms  will  be  defined  as  follows: 


i. 

Corrected  Estimate: 

x  =  x(k  |  k ) 

ii. 

Prediction: 

x  =  x(k  k  - 

in. 

Measurement: 

z 

iv. 

Weight  of  Prediction: 

V. 

Weight  of  Measurement: 

k2 

The  corrected  estimate  can  be  written  as  the  sum  of  the  weighted  prediction  and  the  weighted 
measurement. 


x  =  Kxx  +  K2z 


a) 


In  order  to  proceed,  it  must  be  observed  that: 


X  =  X  +  £ 

z  =  Hx  +  v 


X  =  X  +  £ 


Note  that  s  and  £  are  estimation  errors  with  zero  means,  while  z  is  being  linearized  by  an  H  matrix 
with  additive,  zero  mean  noise.  Hence,  substituting  the  above  equations  into  equation  (1),  we  get: 


X  +  £  —  Kx  (x  +  £  )  +  K2  ( Hx  +  v) 


(2) 


Performing  some  algebraic  manipulation  on  equation  (2)  yields: 
s  =  Kl(x  +  e)  +  K2  (Hx  +  v)  -  x 


=  K^x  +  Kx£  +  K 2 Hx  +  K2v  —  x  (Distributive  Property) 

=  (Krx  +  K2Hx  -  x)  +  K{£  +  K2v  (Commutative  Property) 
=  (Kl  +  K2H  -  I)x  +  Kxs  +  K2v  (Distributive  Property) 


e  =  (Kl+  K2H  -  I)x  +  Kxe  +  K2v 


(3) 
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Recall,  the  means  of  all  the  errors  are  zero  or  in  other  words: 


E[i]  =  E[s]  =  E[v]  =  0 


Therefore,  taking  the  means  of  both  sides  of  equation  (3)  gets: 

E[{KX  +K2H-I)x]  =  0 
E[{KX  +K2H-I)]E[x]  =  0 
(K{  +  K2H  -  I)E[x ]  =  0 

However,  since  E[x\  =  x  the  above  implies: 

Kl  +K2H-I  =  0 


Solving  for  A'  ,  it  is  clear  that: 


Kl  =I-K2H 


(4) 


Let  K 2  =  K  ,  where  K  is  the  Kalman  gain.  Then,  it  is  apparent  that  substituting  equation  (4)  into  equation 
(1)  yields: 


x  =  (/  -  KH)x  +  Kz 
x  =  x  -  KHx  +  Kz 
x  =  x  +  Kz  -  KETx 


x  =  x  +  K(z-Hx) 


(5) 


Hence,  equation  (5)  is  known  as  the  state  correction  equation,  which  corrects  the  prediction  x  by 
weighting  the  measurement  in  a  way  that  minimizes  the  mean  squared  error.  In  order  to  implement 
equation  (5),  it  is  required  that  K  be  expressed  in  known  or  measurable  terms.  In  order  to  do  this  we  must, 
introduce  the  error  covariance  correction: 


P  =  E[££T ] 


In  other  words,  the  error  covariance  correction  is  the  mean  of  the  product  of  the  estimation  error  and  the 
transpose  of  the  estimation  error. 


Recall  that  £  and  £  are  estimation  errors  with  zero  means,  where: 


x  =  x  +  £ 
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X  =  X  +  £ 


Substituting  these  equations  into  (5): 


x  +  £  =  (I-KH)(x  +  £)  +  Kz 
Now,  recall  thatz  =  Hx  +  v ,  yielding: 

x  +  s  =  (I-KH)(x  +  e)  +  K(Hx  +  v) 
x  +  e  =  (I-KH)x  +  (I-KH)s+KHx  +  Kv 
x  +  £  =  x-KHx  +  (I-KH)£+KHx  +  Kv 
£  =  (I-KH)£+Kv 


From  the  definition  of  P  : 

P  =  E{([I  -  KH]£  +  Kv)(£T[I  -  KH]t  +  vtKt)} 


P  =  E{(I-  KH)££T (/-  KH)t  +  (7  - KH)£vtKt  +  Kv£T (/- KH)T  +  KvvT KT } 


P  =  E{(I  - KH)££T (/  -  KH)t }  +  £{(/  - KH)£vtKt)  +  E{Kv£T (7  - KH)T }  +  E{KvvTKT} 
Since  the  estimation  error,  £  ,  is  zero  mean,  and  v  is  zero  mean  non-Gaussian  white  noise: 

E{£VT}  =  E{V£T}  =  0 

Hence,  part  of  the  sum  cancels,  resulting  in: 

P  =  (I-KH)*E{££T}*(I-KH)T  +K*E{vvt}*Kt  (6) 


As  previously  defined: 

P  =  E{££T } 

R  =  E{vvt } 

Therefore,  equation  (6)  becomes: 

P  =  (I  - KH)P(I  -  KH)t  +  KRKt 

Performing  some  algebraic  manipulation  as  follows: 

P  =  (P-  KHP)(I  -  HtKt  )  +  KRKt 
P  =  P-  KHP  -  PHtKt  +  KHPHtKt  +  KRKt 


Now,  taking  the  trace  of  each  side: 

tr(P )  =  tr(P  -  KHP  -PHTKT  +  KHPHT KT  +  KRKT) 


(7) 


(8) 
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This  enables  us  to  combine  some  terms  due  to  the  fact  that: 


tr(KHP)  =  tr(PHTKT) 

Hence,  equation  (8)  is  equivalent  to: 

tr(P)  =  tr(P  -  2 KHP  +  KHPHT KT  +  KRKT)  (9) 

In  order  to  simplify  this  equation  further,  we  can  use  the  fact  that: 


—tr{P)  =  0 
dK 

Therefore,  an  approach  is  to  take  the  derivative  with  respect  to  K  of  both  sides  of  equation  (9): 

—  tr(P)  =  —tr(KHPHTKT  -  2 KHP  +  KRKT) 
dK  dK 

—Hh  =  —HKHPHTKT) — —(2  KHP)  +  —  (KRKt)  (10) 
dK  dK  dK  dK 

For  book  keeping  label  the  parts  of  the  above  equation  as: 

(A)  -Khkhphtkt) 

dK 

(B)  -KtrilKHP) 
dK 


(C)  — HKRKT ) 

dK 


Now,  note  the  following  properties  of - tr  : 

dX 


(i)  XHxyxt)  =  2xy 

dK 


(ii  )Xtr(YXZ)  =  YTZT 
dK 


By  property  (i)  equations  (A)  becomes: 
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-tr{KRKT)  =  2KR 


dK 

Similarly  (C)  becomes: 


,  where  X  is  taken  as  K  and  Y  is  taken  as  R 


- tr{KHPHL  K1  )  =  2 KHPH1  ,  where  X  is  taken  as  K  and  Y  is  taken  as  HPH 1 

dK 


By  property  (ii)  and  the  property  tr(xxT  )  =  tr(xT  x)  equations  (B)  is: 


—tr{2  KHP)  =  —(2  PtKH) 
dK  dK 


- (2 PT KH)  =  2 PHt  ,  where  Y  is  taken  as  PT  and  Z  is  taken  as  H 

dK 

We  now  have  all  the  terms,  transforming  equation  (10)  into: 

—  tr(P )  =  2 KHPHt  -  2 PHt  +  2KR 
dK 


(11) 


Since  tr(P)  is  independent  of  K : 


——tr(P)  =  0 
dK 


Hence,  equation  (11)  results  in: 

2 KHPHt -  2  PHt +  2  KR  =  0 

Equation  (12)  can  be  manipulated  to  finally  yield  K,  the  Kalman  Gain: 


(12) 


2 KHPH1  - 2PH  +  2 KR  =  0 
2 K(HPHt  +R)~  2 PHt  =  0 
K(HPHt  +  R)  =  PHr 


K  =  PHt  (HPHt  +  R)~' 


Now,  let: 


Then  the  Kalman  Gain  is: 


S  =  HPHt  +R 


(13) 


(14) 
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K  =  PHtS  1 


(15) 


Now,  we  will  go  back  to  equation  (7)  in  order  to  express  P  solely  in  terms  of  A’ .  H ,  P ,  and  I  (the 
identity  matrix).  Recall  equation  (7)  is: 

P  =  (I  -  KH)P(I  -  KH)t  +KRKt 

Using  equations  (14)  and  (15),  as  well  as  applying  some  algebraic  manipulation: 

P  =  (P-  KHP)(I  -  KH)t  +  KRKr 

=  P-  KHP  -  PHtKt  +  KHPHtKt  +  KRKt 
=  P-KHP-PHtKt  +  K(HPHtKt  +RKt ) 

=  P-  KHP  -  PHtKt  +  K(HPHt  +  R)Kt 
=  P-  KHP  -  PHtKt  +  KSKt 

=  P-PHTS~XHP-PHT(PHTS~XY  +PHtS~1SS~1HP  (16) 

Make  note  that: 

P  =  PT 

(s~iy  =  s~i 

Then  upon  further  simplification  of  (16),  the  error  covariance  correction  equation  results  as  follows: 

P  =  P  -  PHTS  lHP  - PHTS  lHP  +  KHP 
=  P-  2  PHTS  lHP  +  KHP 
=  P-  2  KHP  +  KHP 
=  P  —  KHP 

P  =  (I-KH)P |  (17) 

Now,  all  that  is  left  is  P ,  which  can  be  derived  from  the  following  set  of  equations: 

(i)  x  -  (f)x 

(ii)  **+i  =  </>xk  +  w 
(ii)  e  =  xk+l  -  x 

The  above  introduces  a  new  matrix  phi,  which  is  the  state  transition  matrix,  and  w,  which  is  a  zero  mean 
white  noise.  Applying  equation  (i)  and  (ii)  in  equation  (iii)  yields: 

s  =  ( (j)x  k  +  w)  -  (j)X 

=  <fi(xk  -  x)  +  w 

=  </>£  +  W 

Recall  that  the  error  covariance  is  the  mean  of  the  estimation  error  and  its  transpose.  Hence,  the  following 
is  the  error  covariance  equation: 
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(18) 


P  =  E{sst } 

=  E  {{(/)£  +  w)((f)£  +  w)T  } 

P  =  </>Pf  +  Q 

Finally,  the  prediction  equation  after  solution  to  differential  equations  assuming  zero-mean,  white  process 
is: 


5c  =  (jk  + 1  P/udt 


+  w 


(19) 


In  summary,  the  Kalman  Filter  equations  are: 


(I) 

x  =  (jk  +  J  +  w 

(II) 

P  =  </>P<f>T  +  Q 

(III) 

K  =  PHt  (HPHt  +  R) 

(IV) 

P  =  (I-  KH)P 

(V) 

x  =  x  +  K(z  -  Hx) 

(VI) 

z  =  Hx  +  v 

2.3  Evaluation 

As  discussed  previously,  the  Kalman  filter  is  an  optimal  filter  for  linear-Gaussian  tracking.  This  technique 
does  not  require  much  space  in  memory  since  only  the  results  of  the  previous  time  step  are  necessary  for 
computation.  The  computation  itself  is  not  complicated,  so  results  are  achieved  rapidly.  It  therefore  finds 
extensive  uses  for  tracking  vehicles  along  roadway  systems,  or  distance-oriented  flight  paths. 

However,  the  Kalman  Filter  is  not  applicable  to  nonlinear  situations  such  as  in  close  combat  or  for  highly 
noisy  measurements  which  might  occur  in  cluttered  environments,  such  as  tracking  ground  transportation 
in  cities.  Therefore  other  techniques  need  to  be  utilized  for  these  situations. 
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3.0  The  Extended  Kalman  Filter 

3.1  Theory 

An  adaptation  to  the  Kalman  Filter  that  addresses  nonlinear  situations  is  the  Extended  Kalman  Filter.  This 
extension  uses  the  same  basic  premise  as  the  KF,  but  it  assumes  the  local  linearity  of  the  track  based  on 
the  definition  of  derivative  in  order  to  prevent  divergence  when  attempting  to  track  highly  nonlinear 
paths.  This  technique  allows  for  more  flexibility  in  its  application.  Although  it  is  not  an  optimal 
estimator,  it  does  approximate  the  state  to  the  first  two  moments  of  its  Taylor  series  expansion. 

3.2  Algorithm 

The  Extended  Kalman  Filter  uses  the  same  equations  as  the  original  Kalman  filter,  replacing  the  state 
transition  matrix  and  the  state-to-measurement  transformation  matrix  with  their  respective  Jacobian 
matrices.  These  Jacobian  matrices  are  as  follows: 

d/i(x)  _  ffl(x) 

dxx  dxn 

:  ;  (39) 

ffmOO  ...  d/m(X) 

dxx  dxn 

,  where  /[/y]  =  Fj-.^or  f[tJ]  =  H[/y] .  In  many  simulations,  the  Jacobian  is  computed  through  the  actual 

derivative  of  the  motion  and  transformation  equations  used  in  the  simulation.  The  definition  of  derivative 
as  the  limit  of  the  following  equation 

/'(*>  = /<Jf  +  Af-/(X>  (40) 

Ax 

as  Ax  goes  to  zero  allows  for  a  numerical  approximation  of  the  Jacobians  by  replacing  a  small  number 
for  Ax  (such  as  1  x  lO-4). 
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3.3  Evaluation 


Unlike,  its  linear  counterpart,  the  Extended  Kalman  Filter,  is  not  an  optimal  nonlinear  filter.  The  EKF  is 
vulnerable  to  divergence  due  to  the  Taylor  Series  approximations  of  its  linearization.  The  EKF  tends  to 
diverge  in  instances  where  there  is  high  initial  uncertainty  in  the  estimate.  This  is  not  to  discount  this 
nonlinear  approximation,  for  in  most  cases  the  EKF  provides  a  sufficient  approximation. 


4.0  The  Particle  Filter 


4.1  Theory 


4.1.1  Exact  Bayesian  Recursion 

The  solution  of  the  filtering  problem  is  to  estimate  the  state  of  a  system  given  all  measurements  up  to 
tim e£ ,  or  in  other  words  to  construct  the  posterior  p.d.f p(xk  \  Zk)  .  Exact  recursive  Bayesian  estimation 

provides  a  solution  in  two  recursive  steps,  prediction  and  update. 

Prediction: 

p(xk  I  zk~ i)  =  J p(xk  I  xk-i)p(xk-i  I  zk-i)dxk-i 


Update: 


p(xk  I  z0  = 


P(Zk  I  )/>(**  I  ZA-l) 

P(Zk  \Zk-X ) 


4.1.2  Sequential  Monte  Carlo  Methods 

Particle  filters  perform  recursive  Bayesian  estimation  directly  on  a  set  of  samples  that  approximate  the 
posterior  density, p{xk  |  Zk)  .  These  techniques,  known  as  Sequential  Monte  Carlo  methods,  differ  from 

exact  recursive  Bayesian  estimation,  in  that  they  do  not  require  the  exact  analytic  solutions  of 
distributions.  An  approximation  of  the  posterior,  achieved  through  the  weighted  sum  of  the  Dirac  Delta 
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function,  is  sufficient  through  which  to  calculate  interesting  statistics,  including  the  minimum  mean 
squared  error  estimate  (MMSE)  of  the  state. 


Monte  Carlo  Integration  &  Importance 

Sampling 

(M.C.  Integration) 

Particle  Filtering  &  Sequential  Monte  Carlo 
Methods  (S.M.C) 

M.C.  Integration  provides  a  method  for 
evaluating  a  multidimensional  integral  by 
sampling  from  a  probability  density  7r(x)  . 

An  estimate  of  the  state  of  a  system  can  be 
attained  from: 

1=  |  g(X)dX  ,  Xe9T 

Elf(xt)\Zk]  =  \f(xtMxt\Zt)dxk 

=  J  /  ( x)n(x)dx 

While,  it  would  be  preferable  to  sample 
directly  from^(x),  this  is  not  always 
possible.  When  sampling  directly  is  not  an 
option,  a  Monte  Carlo  estimate  is  still 
feasible  given  that  an  importance  density  that 
is  similar  to  7r(x)  can  be  easily  sampled  and 
appropriately  weighted. 

In  all  but  the  more  restrictive  cases,  it  is  not 
possible  to  sample  directly  from  the  posterior 
density,  as  a  conceptual  solution  for  the 
propagation  of  this  density  does  not  have  a  closed 
form  solution  or  maybe  be  computationally 
expensive  to  generate  [17].  Therefore,  importance 
sampling  is  applied  and  an  importance  or 
proposal  density  is  chosen: 

/=f fix)  J-  \q(x)dx 

J  q(x) 

The  integral  approximation  achieved  by 
taking  N  independent  samples  from 
q(x)  becomes: 

N 

mx, )  1  Z„)]=j  Ax, )  q(xt  1  z,  )dxt 

q\xk 

Sampling  from  the  proposal  density,  the 
Sequential  Monte  Carlo  estimate  is: 

zi/Wlz,)]=X/(4M4) 

IN  =YJf(xi)w(xi) 

i= 1 

i= 1 

w(xl)  are  the  normalized  importance 
weights: 

q(xl) 

w(xlk  )  are  the  normalized  weights: 

Axi)=p(4lZt) 

9«l  Zk) 

,  w(xl ) 

Ax  )=  ,, 

7=1 

mx‘,a  :ix[) 

7=1 
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4. 1.2.1  Computation  of  Importance  Weights 

The  importance  weights  are  defined  as  a  ratio  of  the  posterior  and  the  importance  density.  In  order  to 
construct  an  implementation  of  the  particle  filter  one  must  simplify  this  ratio  into  quantities  that  are 


known  and  can  be  numerically  computed. 


(I)  Defined  by  Monte  Carlo  Integration: 


w{x[)  = 


p(4  \IA 
|  zk) 


(II)  Simplify  the  joint  posterior  density: 


P(Xk\Zk) 


a.  Applying  Bayes’  Theorem: 

P(Xk\Zk)  = 

b.  Employ  Markov  Chain  Property: 

P(zk  ’  Zk~  i  )p(Xk  ) 


pjZpXQpjXQ 

P(zk) 


P(xk  I Zk)  = 


P(Zk’Zk-l) 


_  Ppk  Iz^-i) 

p(zk  \zk-i ) 


-P(za-  I  Xk,Zk_l')p(xk,Xk_l  |  Zjc-i) 
P(Zk  I  Z k-1 ) 


_  P(.Zk  I  I  Xk_i')p(Xk_l  |  Z^_t) 

P(Zk  I  Z k-\  ) 

c.  The  states  are  propagated  through  a  Markov  process:  p(xk  \  Xk_x )  =  />(x^  |  x^) 
The  observations  are  independent  given  the  states:  p(zk  \  Xk^Zk_l)  =  |  x^) 

The  joint  posterior  reduces  to: 


P(Xk\Zk) 


Ppk  I  **)/>(**  I  *t-l  )/>(**-!  lZA--l) 

/>(**  \zk-x) 


(III)  Combining  (I)  and  (II): 
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w(4 ) 


p(zk  1 4)/K4 1 4-iMVL  I  zi-i) 


4(4  I  Zk)p(zk  \zk~i ) 

(IV)  Let  the  importance  density  factorize  such  that 

<?(4  z*>  =  4(4  1 4-i » zi  )4(4-i  I  ZA'-1 ) 

Then  the  weight  equation  becomes: 

p(zk  1 4  M4 1 4-i )  *  ^(4-i  I  za-i  )  *  1 


w(x'k  )  = 

4(4  I  4-1 »  ZA  )  4(4-1  I  ZA-1 )  P(Zk  I  ZA-1 ) 

(V)  Finally,  it  is  clear  that: 


w( 4)  «  w(4-i) 


pi^i  I4M4 1  4-i  )  I 
4(4  I4-i>za) 


The  mathematical  formulation  above  defines  all  parameters  in  terms  of  probability  densities, 
o  Likelihood  — »  p(zk  \xk) 

o  Transition  Density^  p(xk  \xk-\) 
o  Importance  Density  q(xk  \  xk_x ,  z*  ) 

A  look  at  what  these  densities  mean  in  terms  of  known  statistics  will  reveal  that  these  densities  can  be 
represented  as  integrals  and  are  defined  by  the  system  and  observation  models,  as  well  as  known  noise 
statistics. 

4. 1.2. 1.1  Likelihood 

The  likelihood  characterizes  how  well  the  available  measurement  zk  corresponds  to  the  given 
prediction  xk .  “The  conditional  PDF  of  zk  given  xk,  p(zk  \xk),  is  defined  by  the  measurement 
model  and  the  known  statistics  of  observation  noise,  vk  ”  [10]. 

p(zk  I  xk )  =  J  5(zk  -  K  (xk ,  vk  ))p(vk  )dvk 
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The  equation  hk  (xk ,  vk )  is  known  as  the  observation  equation.  This  function,  which  is  a  known  model, 
relates  the  measurements  to  the  state  vector. 


Derivation  of  Likelihood  Equation 

(I)  Integration  allows  the  likelihood  to  be  written  as: 

p(zk  xk )  =  j \p(zk’vk  I  xk)dvk 

(II)  Applying  Bayes’  Theorem: 

•p(.zk>vk’xk) 


p(zk  I  xk)= 


p(xk ) 


-d\u 


(III)  Implementing  the  Markov  Chain  Rule  produces  the  following  manipulation: 

p(zk  i  xk ) = j p(zk  i  vk’xk )p(vk  i  xk )dvk 

(IV)  Assume  vk  is  independent  of  xk ,  that  is  the  observation  noise  is  independent  from  past  and 
current  states,  then  p(vk  \xk)  =  piyk)  and  the  likelihood  becomes  equivalent  to: 

p(zk\xk)  =\p(zk\vk’xk)p(vk)dvk 

(V)  Manipulating  as  to  incorporate  the  Dirac  Delta  function: 

Note  that: 


P(zk  I  vk’xk)  = 


' 1  when  hk{vk,xk)  =  zk 
0  when  hk(vk,xk)^z 


Hence,  it  follows  that: 


lim 5(zk  ~ K (vk >xk))  =  p(zk  I  vk,xk) 


Applying  the  above  notion  to  (IV)  gives  the  following  result: 


p(zk  I  xk )  =  J  5izk  -  K  (xk  >  vk  ))p(yk  )dvk 


The  computation  of  the  likelihood  is  dependent  upon  the  probability  of  observation 
noise,  as  indicated  by  the  derived  equation. 
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There  may  not  generally  be  only  a  concern  with  utilizing  information  from  a  single  sensor.  The  extension 
to  a  multiple  sensor  problem  is  not  difficult.  Under  the  assumption  of  N  independent  sensors,  centralized 
measurement  fusion  in  the  PF  is  simply  computed  by: 

sN 

n  p(zjk  xk ) 

j= i 

4. 1.2. 1.2  Transition  Prior 

The  transition  prior  describes  the  evolution  of  the  state.  Where  the  likelihood  is  defined  by  the 
observation  equation  and  known  measurement  noise,  the  transition  prior  p{xk  \  xk_x)  is  defined  by  the 

known  system  equation  fk_x  ( xk_x ,  wk_x )  and  process  noise  wk  .  A  mathematical  representation  of  the 
transition  prior,  which  can  be  derived  in  the  same  manner  as  the  likelihood,  takes  on  the  form: 


p(xk  i  xk- 1 ) = J  s(xk  -  fk~i  (v-i > wk~  i  ))p(wk~  i  yiwk~< 


4. 1.2. 1.3  Importance  Density 

The  proposal  density  q(xk  \  xk_x,yk)  is  a  design  choice  involved  in  the  implementation  of  a  particle 

filter.  In  choosing  a  proposal  density  one  would  want  a  density  as  similar  to  the  posterior  as  possible. 
However,  the  Monte  Carlo  Method  only  requires  that  the  proposal  and  posterior  have  the  same  support. 
That  is,  particles  that  are  associated  with  a  probability  greater  than  zero  in  the  posterior  should  also  be 
greater  than  zero  in  the  proposal. 

4. 1.2. 1.3.1  Transition  Prior  Proposal 

The  most  commonly  used  proposal  distribution  is  the  transition  prior  p{xk  Selection  of  the 

transition  prior  as  the  proposal  density  transforms  the  generic  Sequential  Importance  Sampling  algorithm 
into  the  algorithm  known  as  the  Bootstrap  Filter.  The  prior  is  a  popular  choice,  for  it  is  easy  to  implement, 
causing  the  weight  computation  to  reduce  to  a  computation  of  the  likelihood. 


wk  =  Wk-lPpk  I  Xk) 

The  decision  to  choose  the  prior  as  the  sampling  proposal  disregards  the  fact  that  the  proposal  should  take 
into  account  the  most  recent  observation.  As  a  result,  the  filter  may  potentially  breakdown  when  the  prior 
and  likelihood  do  not  overlap.  Essentially,  the  likelihood  weights  improbable  particles  more  highly  than 
probable  particles. 
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Prior 


Likelihood 


*  •  •  ••••••  t  i  •  •  •  i 

Above  image  is  from  http://cslu.cse.ogi.edu/publications/ps/UPF  CSLU  talk.pdf .  This  image  illustrates 
the  likelihood  and  prior  disagreement  that  can  lead  to  sample  degeneracy. 

4. 1.2. 1.3. 2  Alternative  Proposals 

Alternatives  to  the  transition  prior  proposal  include  using  other  fusion  algorithms  as  the  proposal  density, 
including  the  EKF  and  Unscented  Kalman  Filter  (UKF).  These  approaches  may  be  favored,  as  they  take 
into  account  the  current  measurement.  However,  using  these  proposals  is  both  computationally  more 
complex  and  as  a  result  requires  longer  processing  time. 


4. 1.2. 2  Resampling 

Resampling  is  the  final  step  of  an  iteration  of  the  particle  filter.  Without  the  resampling  step,  the 
algorithm  so  far  derived  is  known  as  Sequential  Importance  Sampling  (SIS).  The  resampling  step 
becomes  a  necessity  due  to  the  fact  that  one  must  choose  an  importance  density  and  not  sample  directly 
from  the  posterior.  As  a  result  of  the  inability  to  sample  directly  from  the  posterior,  the  variance  of  the 
weights  increases  over  time,  resulting  in  a  phenomenon  known  as  the  Degeneracy  Phenomenon. 
Essentially,  this  phenomenon  causes  the  number  of  effective  particles  to  decrease  over  time  and  may  after 
many  iterations  result  in  only  one  particle  with  non-negligible  weight.  Hence,  to  prevent  this  degeneration 
resampling  reselects  particles  with  high  importance  weights.  The  following  figure  illustrates  the 
resampling  method  for  the  reselection  of  particles  with  higher  importance  weights. 


4. 1.2. 2.1  Resampling  Schedules 

(I)  Deterministic  resampling  describes  a  schedule  with  fixed,  interval  times  ( tx ,  t2 ). 
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(II)  Dynamic  scheduling  is  done  when  the  effective  number  of  particles  is  less  than  a  predetermined 
threshold  number. 

4. 1.2. 2.1  Resampling  Schemes 

4. 1.2. 2. 1.1  Systematic  Resampling 

Systematic  Resampling  is  the  computationally  easiest  resampling  scheme.  This  ease  comes  from  the  fact 
that  it  uses  the  simple  uniform  distribution.  The  following  equation  is  implemented  when  performing 
systematic  resampling: 


1)  Ul  =  - — -  +  U ,  where  U  is  a  single  random  drawn  from  ((0,1/n])  and  /  = 


Th qU1  is  the  ith  uniform  sample.  That  is,  the  above  equation  randomly  describes  n  uniform  samples  on 
the  interval  (0,1].  This  sample  is  created  by  first  drawing  a  random  number  on  the  interval  (0,1/n].  The 

i  —  1 

remaining  uniform  samples  are  then  generated  by  adding  -  to  the  randomly  sampled  U.  Now,  re- 

n 


index  the  jth  particles  based  on  the  n  ith  uniforms. 


Even  though  the  systematic  resampling  scheme  is  most  popular  due  to  the  easy  computation  of  a  uniform 
distribution,  other  methods  of  resampling  have  been  developed.  The  main  differences  between  these 
resampling  schemes  are  the  distributions  used  to  re-index  the  sample  set. 


4. 1.2. 2. 1.2  Multinomial  Resampling 

In  multinomial  resampling,  the  duplication  count  of  the  jth  particle  is  determined  by  a  multinomial 
distribution.  That  is,  Nx , . . . , N m  are  distributed  according  to  the  multinomial  (n;  wx , . . . ,  wm  ). 


Process: 

1)  Draw  n  independent  uniforms  {Ul  },  1  <  i  <  n  on  the  interval  (0,1]. 

ri'  rV  r  D inv  (U* ) 

2)  Set  1 1  —  D  (U  1  )  and  9  —  —  ,  where  D™v  is  the  inverse  of 

the  cumulative  distribution  function  associated  with  normalized  weight,  that  is, 

D  7  ( u  )  =  i  for  u  e  (]T  wJ ,  2  ] . 


7  =  1 


7=1 
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4. 1.2. 2. 1.3  Residual  Resampling  (Remainder  Resampling) 

1)  Residual  resampling  is  determined  by  the  following  equation: 

Nj  =  [nwj  \+Nj 

NJ :  total  number  of  times  the  jth  particle  is  duplicated 

Nj  :  generated  according  to  a  multinomial  distribution  Mult(n  -  R;wl  , . . . ,  wn  ) 

R  :  the  sum  of  the  integer  part  of  the  product  of  the  total  number  of  particles  and  the  weight  of 
the 

jth  particle  (XM) 

— 7  nwj  -  nwj 

wJ  :  - - - -a  =  1, ...,n 

n-R 

\nwJ  J:  integer  part  of  the  product  of  the  total  number  of  particles  and  the  weight  of  the  jth 
particle. 


4.2  Generic  Particle  Filter  (Pseudo  code) 

•  Initialization: 

-  Draw  N  particles  from  a  known  initial  distribution  p(xx ) 

•  Update: 

-  Upon  receiving  a  measurement,  evaluate  the  importance  weights  according  to: 

,  p(y\  I  A)p(x\) 

—  - : - 

tfOi'IJi) 

-  Normalize  importance  weights: 

p(y 1 1 4)p(xl) 

ypjyi  IV MV) 
j=\  q(x(  I  j;i  ) 

-  Calculate  Number  of  Effective  Particles 

-  If  the  number  of  effective  particles  falls  below  a  threshold  (IF  Neff  { Nthr ) 

Resample  using  a  chosen  Algorithm 

•  Prediction: 

FOR  i  =  l:N 
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-  Draw  xlk  ~  q(x[  \xlk_l,yk)by  passing  particles  obtained  after  resampling  through 
proposal  density  (Note:  In  the  bootstrap  filter,  one  would  pass  the  particles  through  the 
system  equation  to  draw  these  x\  particles.) 


•  Update: 

-  Upon  receiving  a  measurement,  evaluate  the  importance  weights  according  to: 

p(yk\xk)p(.xk\xk-i) 

Wk  ~  - : - : - 

^(VK-n.y*) 

-  Normalize  importance  weights: 

pjyjc  |4M4  I4-i) 

-/  qtfk  1 4-i >yk) 

Wk  =  — - : - : - : - 

yPin  I  xJk)P(xJk  K-i) 

y'=l  UwJt) 

-  Calculate  the  number  of  effective  particles 

-  If  the  number  of  effective  particles  falls  below  a  threshold  (IF  Neff  (Nthr ) 

Resample  using  a  chosen  Algorithm 


5.0  Experiments 

5.1  Maneuvering  Targets 

The  algorithms  described  thus  far  are  implementations  using  a  single  dynamic  model.  For  tracking 
applications,  a  single  model  algorithm  may  not  fully  capture  the  true  motion  of  a  target.  While,  targets 
may  move  with  a  constant  velocity  for  a  period  of  time,  a  target  may  at  a  different  time  execute  a  turn,  or 
in  other  words  maneuver.  In  order  to  account  for  the  varied  motion,  algorithms  that  implement  multiple 
modes  that  are  modeled  as  a  Jump  Markov  System  (JMS)  are  used. 

5.2  System  Models 

The  motion  of  the  target  is  described  at  a  given  time  instance  by  one  of  three  dynamic  models.  The  three 
models  of  target  motion  used  are:  (1)  Constant  Velocity  Model  (CV)  (2)  Clockwise  Coordinated  Turn  Model 
(CT)  (3)  Counterclockwise  CT  Model.  These  three  models  are  used  in  a  switching  manner  to  simulate  a 
maneuvering  target,  where  a  regime  variable 

rk  e  { 1,  2,  3  1 1  =  CV,  2  =  Clockwise  CT ,3  =  Counterclockwise  CT} 
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is  used  to  account  for  the  current  mode  at  discrete  time  index  k.  The  switching  between  models  follows 
transitions  given  by  a  Markov  chain  with  transition  probabilities  Ktj  =Pr{^+1  =j\rk  =/},(/,/ e{l, 2,3}),  such 

that  Ky  >  0,  —  1 .  The  evolution  of  the  states  xk  can  then  be  described  as 

j 

=f(rt\xk_l)  +  Gwk 

where  the  state  vector  x  k  is  defined  as  x ,  =  [.v  y  x  y\,wk  is  the  Gaussian  process  noise  with 
covariance  matrix  Q  =  q*  I ,  where  /  is  an  identity  matrix  and  q  is  a  scalar  ,  and 


H  0 
2 

T 2 

G=  0  — 

2 

T  0 

0  T 

with  T  defined  as  the  sampling  interval.  The  discrete  time  instances  where  the  state  evolves  according  to 
the  CV  model  with  regime  r,  =  1  can  be  gotten  by  replacing  f{l>:  1  )  with  the  transition  matrix: 


Fm(*k)  = 


1  0  T 
0  1  0 
0  0  1 
0  0  0 


0 

T 

0 

1 


For  evolution  of  the  state  with  a  constant  turn  rate  the  transition  matrices  are 


1  0 


Fu\xk)  = 


0  1 


0  0 
0  0 


sin(Oy  *  T) 

op 

(l-cos(OP  *T) 

op 

cos (Okj)  *  T) 
sin(0^)  *  T) 


-  (1  -  cos(0^/)  *  T )) 

op 

sin(C)^)  *  T) 

op 

-  sin(Oy,)  *  T ) 
cos iPJk  *  T) 


,7=2,3 


with  clockwise  and  counterclockwise  turn  rates  given  by 


oP  = 


-fx 


v2  +y2 
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of  =  - 

where  am  is  a  constant  maneuver  acceleration  parameter. 

5.2.1  Interacting  Multiple  Model  EKF  (IMM-EKF) 

The  multiple  dynamic  models  are  implemented  in  the  IMM-EKF  as  a  process  of  interaction,  filtering,  and 
combination  of  outputs.  The  interaction  begins  each  iteration  of  the  algorithm  by  mixing  the  outputs  from 
the  previous  step  based  upon  the  previous  mode  probabilities  .  Next,  the  inputs,  derived  from  a  single 

EKF  with  associated  dynamic  model,  are  filtered,  resulting  in  both  updated  model  probabilities  and  new 
inputs  for  the  next  step.  Finally,  the  current  model  probabilities  are  used  to  compute  a  weighted  average 
of  the  outputs  from  the  individual  filters,  resulting  in  a  combination  and  an  estimate  for  discrete  time 
step  k . 

5.2.2  Multiple  Model  PF  (MM-PF) 

The  alterations  in  the  MM-PF  compared  to  the  original  algorithm  are  not  drastic  and  do  not  significantly 
increase  the  computational  complexity  of  the  regular  SIR  algorithm.  The  added  computations  involve  the 
determination  of  the  regime  variable  for  each  of  the  n  particles  at  discrete  time  step  k .  The  mode  is 
computed  based  upon:  if  rk_x=i  and  un~U[ 0,1],  then  {m  e  {1,2,3}  |  rk  =  m}  such  that: 

m— 1  m 

^^71-  <un<  ^71-  .  The  particles  are  then  propagated  through  the  filter  and  predictions  are 

7=1  7=1 

determined  that  take  into  account  the  regime  variable  selected. 

5.3  Measurement  Models 

The  implemented  simulation  employs  the  ability  to  select  sensors  of  four  types  that  provide  the  following 
types  of  measurements:  1)  Bearings  Only  2)  Range  Only  3)  Range  and  Bearings  4)  Range,  Elevation,  and 
Azimuth.  The  measurement  model  described  as 

z{  =  hu\xk)  +  vk 

where  vk  is  the  measurement  noise,  can  employ  any  of  j  sensors  with  /  =  1,...,  4 measurement  models. 
The  measurements  models  for  the  sensors  can  take  any  of  the  following  forms: 
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where  ( xSj ,  ySj  ,zSj)  is  the  position  of  sensor  j  . 

5.4  Measures  of  Performance 

Evaluations  of  the  PF  and  EKF  algorithms  were  primarily  based  upon  average  mean  squared  error  (MSE) 
and  variance  (V)  over  separate  Monte  Carlo  Runs.  The  mean  squared  error  and  variance  computations 
are: 

#  Time  Steps  #Runs  . - 

E  E  ~xU2 

MSE  ,  =  — — - J— - 

(#  Runs  )(#  TimeSteps  ) 


WTimeSteps  #Runs  i - ; - 

E  Ev(x‘~-<„)2 

^=1  i-l _ 

(  #Runs  )(  #  TimeSteps  ) 
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5.5  Simulations 


Comparison  of  the  EKF  vs.  PF  depends  on  several  factors  including:  sensor  position,  “q”,  and  number  of 
particles  used  in  the  particle  filter. 

Sensor  position  can  degrade  overall  tracking  performance  for  both  EKF  and  PF.  Overall  tracking 
performance  refers  to  RMSE  and  variance.  A  high  tracking  performance  is  described  as  having  low 
RMSE  and  low  variance.  Performance  is  degraded  when  the  sensors  are  not  lined  up  normal  to  the  target 
track.  When  the  sensors  are  lined  up  in  parallel  with  the  direction  of  motion  of  the  target,  performance  is 
minimized. 

In  this  simulation,  overall  tracking  performance  was  at  a  minimum  when:  1)  the  two  sensors  were 
positioned  side  by  side  on  all  four  sides  of  the  target’s  track;  2)  the  two  sensors  were  positioned  on 
opposite  sides  of  the  target  track  parallel  with  the  target’s  motion.  Alternatively,  tracking  performance 
was  maximized  when  either:  1)  one  or  more  of  the  sensors  were  positioned  close  to  the  target  track  so  at 
no  time  during  the  simulation  was  the  sensor  in  parallel  with  the  target’s  motion;  2)  the  two  sensors  were 
positioned  at  right  angles  on  the  outside  of  the  target  tack  ensuring  at  least  one  sensor  was  normal  to  the 
target  at  all  times. 

In  the  two  cases  mentioned  previously  where  performance  is  minimized,  the  EKF  outperformed  the  PF  in 
nearly  every  situation.  Despite  the  EKF’s  performance  over  the  PF  in  this  sensor  setup,  it  is  important  to 
remember  the  RMSE  and  the  variance  were  still  very  large. 

When  the  sensors  are  positioned  so  that  performance  is  maximized  the  results  weren’t  as  one-sided 
toward  the  EKF  as  when  the  sensors  were  positioned  poorly.  A  noticeable  trend  emerged  from  this  sensor 

setup.  The  PF  outperformed  the  EKF  whenever  “q”  was  set  to  1*6  x  10  no  matter  how  many  particles 
were  used  (50,  100,  500). 

When  the  PF  was  set  to  50  particles  in  this  performance  maximization  sensor  setup,  the  EKF  showed 
better  performance  for  the  “q”  values  of  0.5,  1,  and  5.  As  the  number  of  particles  increased,  the 
performance  of  the  PF  increased  (with  the  exception  that  simulation  time  increased  dramatically  as  a 
result  from  more  computational  complexity).  When  the  number  of  particles  was  increased  to  100  the 

EKF  outperformed  or  did  just  as  well  as  the  PF  with  the  exception  of  when  “q”  was  set  to  1  ^  x  10  as 
mentioned  above.  Once  the  number  of  particles  was  set  to  500  the  PF  either  outperformed  or  did  just  as 
well  as  the  EKF. 
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2  Positron  (km) 


Examples  of  Poor  Sensor  Alignment 
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2  Position  (km)  Z  Position  (km)  Z  Position  (km) 


Examples  of  Proper  Sensor  Alignment 


Figure  3.  At  least  one  sensor  positioned  close  to  the  track  path. 


Figure  5.  Sensors  positioned  at  right  angles  centered  on  the  target  track. 
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Typical  Results  with  Poor  Sensor  Alignment 


Sensor  Positions  and  True  Trajectory 


Figure  6.  Average  tracks  of  100  Monte  Carlo  simulations  displaying  typical  results  with  the  two  sensors 
positioned  on  opposite  sides  of  the  target  track. 


Average  Tracks  of  PF  &  EKF  v  True 


Figure  7.  Average  tracks  of  100  Monte  Carlo  simulations  displaying  typical  results  with  the  two  sensors 
positioned  side  by  side. 
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Effects  of  “q”  and  the  Number  of  Particles  with  Properly  Aligned  Sensors 


O 

CL 


Figure  8.  Average  tracks  of  100  Monte  Carlo  simulations  with  q  =  1.6x10  6  and  the  number  of  particles  set 


to  50. 


X  Position  (km) 

Figure  9.  Average  tracks  of  100  Monte  Carlo  simulations  with  q  =  5  and  the  number  of  particles  set  to  50. 
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Y  Position  (km) 


Average  Tracks  of  PF  &  EKF  v  True 


Figure  10.  Average  tracks  of  100  Monte  Carlo  simulations  with  q  =  1.6x10  6  and  the  number  of  particles 
set  to  100. 


Average  Tracks  of  PF  &  EKF  v  True 


X  Position  (km) 

Figure  11.  Average  tracks  of  100  Monte  Carlo  simulations  with  q  =  5  and  the  number  of  particles  set  to  100. 
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Average  Tracks  of  PF  &  EKF  v  True 


Figure  12.  Average  tracks  of  100  Monte  Carlo  simulations  with  q  =  1.6x10  6  and  the  number  of  particles 
set  to  500. 


Average  Tracks  of  PF  &  EKF  v  True 


Figure  13.  Average  tracks  of  100  Monte  Carlo  simulations  with  q  =  5  and  the  number  of  particles  set  to  500. 
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6.0  Chief  Scientist  Briefing  Results 


On  August  30,  2007,  a  briefing  was  given  by  Mark  Alford  to  the  Chief  Scientist  of  AFRL/IF,  Dr.  John  S. 
Bay.  The  following  is  a  summary  of  that  briefing. 

The  overall  in-house  program  is  being  referred  to  as  “Fusion  Techniques  and  Non-Linear  Filtering 
(FTNLF).”  This  overall  program  encompasses  past  and  future  in-house  efforts  along  with  component 
contractual  efforts.  The  outline  of  the  briefing  consisted  of  going  over  the  purpose,  particle  filtering,  and 
progress,  followed  by  a  discussion  of  the  major  in-house  efforts  and  then  the  component  efforts,  followed 
by  conclusions  and  future  directions. 

The  primary  purpose  of  the  FTLNF  In-House  Program  is  to  investigate  and  develop  promising  innovative 
technologies  that  hold  promise  for  improvements  in  Air  Force  target  tracking  and  multi-sensor  fusion 
capabilities.  Current  permanent  staffing  is  shown  in  the  following  organizational  chart. 


Figure  14.  In-House  FTLNF  program  staffing 


Also,  there  are  typically  two  summer  students  hired  each  year  to  help  support  this  program.  In  the 
summer  of  2007,  Christopher  Poore  and  Becky  Bailey  worked  on  nonlinear  filtering  techniques,  and  their 
reports  are  included  in  the  Appendices  of  this  report.  There  is  also  University/Contractual  support 
through  Dr.  Pramod  Varshney,  Syracuse  University,  Black  River  Systems  Company,  and  Numerica 
Corporation. 


Figure  15  shows  an  overview  of  the  fusion  scenario. 
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Figure  15.  Overview  of  the  Fusion  Scenario 

Multiple  assets  are  shown  conducting  surveillance  over  the  ground,  air,  space  and  sea.  Some  have 
overlapping  coverage,  some  do  not.  This  shows  the  idea  of  sharing  information  between  platforms  and 
with  ground  and  space  based  systems.  These  platforms  and  systems  are  monitoring  events  in  both  war 
and  peace  time  environments.  Examples  of  events  are  the  fires  on  the  ground.  Terrain  features  are  also 
important  in  this  scenario.  Both  centralized  and  decentralized  fusion  centers  are  being  employed. 

The  problem  statement  is  that  Nonlinear  Non-Gaussian  Processes  (NNGP)  present  a  major  challenge  in 
all  types  of  military  problems.  This  is  because  the  real-world  is  nonlinear  and  non-Gaussian,  despite 
assumptions  made  in  most  conventional  fusion  algorithms.  The  FTLNF  In-House  program  is  addressing 
this  problem  by  researching  and  developing  tracking  filters  that  do  not  presume  a  linear  Gaussian  world. 
Perhaps  the  most  famous  of  these  is  the  particle  filter.  The  theory  behind  particle  filtering  has  been 
around  for  a  long  time,  namely  in  the  form  of  Monte  Carlo  Markov  Chain  based  approaches,  however,  it 
resurfaced  with  the  name  particle  filtering  in  1993  with  Neil  Gordon’s  thesis  referenced  in  [10]. 
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Figure  16  depicts  an  urban  clutter  picture  which  shows  how  nonlinear  things  can  be.  Imagine  a  target 
moving  through  this  environment.  Clearly  this  is  a  nonlinear  environment.  With  a  little  imagination,  you 
can  almost  tell  that  it  is  non-Gaussian  as  well.  Consider  a  camera,  radar,  infrared  sensor,  unmanned  aerial 
vehicle  (UAV),  or  other  sensing  device  overlooking  this  scenario.  The  question  is,  how  do  we  discern 
critical  situations  in  this  sort  of  environment.  This  is  exactly  where  it  is  expected  that  new  nonlinear 
filtering  techniques  will  help. 


Figure  16.  Urban  Clutter  Nonlinear  Environment 


The  particle  filter  is  the  focus  of  the  FTNLF  research.  Particle  Filtering  holds  great  promise  to  provide 
enhanced  capability  for  non-linear,  non-Gaussian  tracking  conditions.  Considerable  work  is  required  to 
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verify  the  situations  where  the  Particle  Filter  will  provide  improved  performance.  Developing  tools  and 
techniques  to  investigate  Particle  Filter  performance  relative  to  other  filtering  techniques  is  a  major 
emphasis  of  the  effort  under  the  FTNLF  In-House  program.  Other  non-linear  filters  are  also  being 
investigated  for  comparison  purposes  to  determine  the  conditions  under  which  they  represent  a  better 
choice  than  the  particle  filter. 
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Figure  17.  Particle  Filter:  A  Block  Diagram 

Figure  17  shows  a  block  diagram  of  the  particle  filter.  The  input  is  the  probability  density  function  (pdf) 
of  the  target  state  x  consisting  of  position,  velocity,  acceleration,  jerk  and  other  attributes.  The  pdf  goes 
through  the  state  update  equation  to  propagate  the  initial  state  (prediction  step)  from  previous  time  step  k- 
1  up  to  time  k.  A  measurement  comes  in  consisting  of,  for  example,  range-bearing  or  bearings  only  from 
a  radar,  infrared,  or  Electronic  Support  Measures  (ESM)  sensor.  This  measurement  is  passed  through  the 
measurement  equation  to  form  a  likelihood  that  measurement  y  was  generated  by  state  x.  The  weighted, 
updated  pdf  is  then  resampled.  When  this  measurement  updated  density  is  formed,  the  mean  is  calculated 
to  create  the  track.  The  whole  process  is  then  repeated  for  the  next  time  step  and  measurement  update. 

The  FTNLF  In-House  program  progress  in  Particle  Filter  investigation,  development,  and  verification  has 
focused  on  the  development  and  utilization  of  a  simulation  capability  including: 

1 .  Development  of  a  MATLAB  simulation  capability  in  coordination  with  Syracuse  University 

2.  Preliminary  investigation  of  Particle  Filtering  capabilities 
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Develop  Particle  Filter,  - 
compare  to  Kalman  Filter, 
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scenarios  and  evaluation. 


-Scenario  development. 
-Manually  swap  filtering 
modules  and  perform 
test  and  evaluation  in 
more  complex  and 
dynamic  situations. 
Minigrant  #2 
Image  Based 
Particle  Filter _ 

Numerica  -  Particle 
Filter  to  track  objects 


2)  Develop  and  automate 
Interacting  Multiple  Models 
for  filter  swapping  based 
upon  situation 

3)  Enhance/mature  Image 
Based  Particle  Filter 
from  Minigrant  #2 


in  video 


Figure  18.  FTNLF  In-House  Program  Timeline 


Figure  18  shows  a  timeline  of  the  in-house  programs  comprising  FTNLF.  The  first  program  was 
Complementary  Advanced  Fusion  Exploration  (CAFE)  which  was  a  study  to  identify  promising  types  of 
algorithms  that  could  aid  in  problems  of  advanced  fusion  and  tracking.  It  was  found  that  Particle  Filtering 
was  an  area  that  was  not  being  worked  within  the  Information  Directorate  at  Rome,  and  was  chosen  as  an 
area  for  detailed  examination.  A  final  report  was  published  for  CAFE  and  a  small  examination  of  particle 
filtering  was  done.  It  showed  that  the  particle  filter  outperformed  the  extended  Kalman  filter  for 
nonlinear  trajectories  and  Rayleigh  Noise.  In  fact,  it  was  found  that  the  Unscented  Particle  Filter  (UPF) 
outperformed  the  Unscented  Kalman  Filter  (UKF)  and  the  Particle  Filter. 


The  follow-on  to  CAFE  was  Development  and  Evaluation  of  Fusion  Techniques  (DEFT)  whose  objective 
was  to  understand  the  Particle  Filter  in  more  depth.  The  result  was  the  development  of  a  particle  filter  to 
compare  to  Extended  Kalman  filters  with  some  preliminary  scenarios  and  evaluation.  The  Multi-INT 
Particle  Filter  Minigrant  was  an  outgrowth  of  DEFT.  DEFT  developed  a  two-dimensional  particle  filter 
for  in-house  testing  and  analysis.  Other  programs  initiated  under  DEFT  were  the  In-House  Particle  Filter 
Analysis  and  Testing  (IHPFAT)  fallout  money  program  and  another  Minigrant  to  do  Image  Based 
Particle  Filtering.  IHPFAT  has  the  objective  of  combining  the  kinematic  and  particle  filter  with  the 
Baseline  Road  Assisted  Tracker  developed  by  Black  River  Systems  Company.  Additionally,  Black  River 
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was  tasked  with  helping  AFRL  with  appropriate  scenario  development.  Once  these  software  tools  are 
integrated,  the  idea  was  to  manually  swap  filtering  modules  and  perform  test  and  evaluation  in  more 
complex  and  dynamic  situations.  The  Image  Based  PF  Minigrant  is  just  getting  under  way,  so  not  much 
progress  has  been  made,  but  the  idea  is  to  use  a  particle  filter  to  track  objects  in  video. 

Tracking  Evasive  Nonlinear  Targets  is  a  new  two  year  In-House  program  to  perform  analysis  and 
enhancements  based  on  lessons  learned  under  previous  efforts.  The  ultimate  goal  is  to  develop  and 
automate  an  Adaptive  Nonlinear  Tracking  System  (ANTS)  for  filter  swapping/parallelization  based  upon 
situation.  It  is  also  a  goal  to  enhance/mature  Image  Based  Particle  Filtering  from  the  second  Minigrant. 


Under  DEFT  there  were  some  accomplishments  that  lead  to  the  conclusion  that  particle  filtering  has 
promise.  DEFT  developed  a  two-dimensional  particle  filter  for  in-house  analysis  and  testing.  Under 
DEFT,  the  team  performed  Monte  Carlo  simulation  runs  to  test  the  performance  of  the  Particle  Filter  (PF) 
as  compared  to  the  EKF  for  multiple  bearings  only  sensors  (ESM  sensors).  A  comparison  was  made  of 
tracking  error  variance  based  on  Root  Mean  Square  (RMS)  position  error.  Sample  results  (where  the  PF 
outperformed  the  EKF)  are  shown  in  Figures  19  and  20. 


Experiment  007 


Experiment  number  7  in  series  of  15. 

Number  of  particles:  50 

Number  of  runs:  100 

Number  of  time  steps:  60 

Active  Sensors  Type  x  pos  y  pos  glint 

probability 

4  1  30  -50  0.0 

5  1  0  -100  0.0 

Number  of  active  sensors:  2 
Maneuverability  variable  q:  1.600e-005 

PF  lost  track  0  times  out  of  1 00,  or  0.00%  of  the 
time. 

EKF  lost  track  28  times  out  of  1 00,  or  28.00%  of  the 
time. 

AvgTime/Time_Step  Avg  RMSE  Variance 

PF  0.0417  1.795  3.000 

EKF  0.0042  3.449  47.631 

Processing  took  9.91  times  longer  for  PF  than  for 
EKF. 


Figure  19.  DEFT  Simulation  Conditions:  2  ESM  Sensors 
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Exp007:  EKF  -  All  Tracks 


Exp007:  Position  Variance  PF  v  EKF 


8-4o  -40  -30  -20  -10  0  10 

X  Position  (km) 


Exp007:  PF  -  All  Tracks 


much  less  for  PF  than  EKF 


Figure  20.  Results:  variance  is 


Figure  7  shows  that  the  tracking  error  variance  is  much  better  for  the  PF  than  the  EKF,  requiring  lOx  as 
much  processing  time.  The  result  is  for  50  particles,  100  Monte  Carlo  runs,  two  bearings  only  sensors. 
The  EKF  lost  track  28%  of  the  time  whereas  the  PF  did  not  lose  track  at  all.  The  variance  and  RMSE 
versus  average  computation  time  per  time  step  are  compared  in  Figure  6.  Use  of  High  Performance 
Computers  (HPCs)  is  being  worked  with  IFTC. 


The  follow  on  effort  to  DEFT  is  TENT.  TENT  will  develop  tools  and  techniques  such  as  automated  filter 
selection  and  improved  tracking  for  nonlinear:  target  motion,  platform  motion,  and  sensor  dynamics. 
TENT  will  investigate  the  applicability  of  image  based  tracking,  enhance  existing  algorithms,  and 
leverage  contributing  technologies  such  as  Multi-INT  Particle  Filtering. 
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Tracking  Evasive  Nonlinear  Targets 

(TENT) 


Tracking  Targets  in  Non-linear  Motion 


Description: 

This  in-house  effort  will  develop  tools  and 
techniques  to  automate  filter  selection  to  improve 
tracking  by  reducing  susceptibility  to  non- 
linearities  in  target  motion,  platform  motion,  and 
sensor  dynamics.  It  will  also  investigate  the 
applicability  of  image  based  particle  filter  tracking 
as  well  as  enhancing  existing  algorithms. 


Literature  Search 

Investigation  (Conceptualization) 

Software  Development 

Analysis  (Proof  of  Concept) 

Laboratory/Metrics  definition 

Experimentation 

Burdened  Funding 
(K) 

Benefits  to  the  War  Fighter: 

•  Provides  new  capabilities  for  tracking  targets 
moving  in  nonlinear,  non-Gaussian  manner 

•  Promotes  longer  tracks  by  having  less  target  loss 
during  maneuvers  (fewer  dropped  tracks) 

•  Dismount  Tracking 
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Technologies:  Multi-INT  Particle  Filtering 


Figure  21.  TENT  Quad  Chart 


A  quad  chart  for  TENT  is  shown  in  Figure  21.  The  TENT  objectives  are  Nonlinear  Development  and 
Enhancements,  Tracker  that  employs  a  combination  of  nonlinear  filters  and  enhance  and  mature  image 
based  tracking.  The  overarching  theme  is  to  increase  track  lifetimes  by  developing  fusion  techniques 
resilient  to  evasive  maneuvers. 
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TENT  Objectives 


*  Nonlinear  Filtering  Development  and  Enhancements 

—  Analyze  effects  of  update  rates  and  realistic  scenarios  that 
present  varying  degrees  of  nonlinearities 

®  Tracker  that  Employs  a  Combination  of  Nonlinear  Filters 

—  Enhance  nonlinear  filtering  by  adaptively  switching  or 

parallelizing  between  algorithms  based  upon  characteristics, 
performance  and  constraints  of  a  situation 

#  Enhance  and  Mature  Image  Based  Tracking 

—  Integrate  image  particle  filtering  with  existing  Multi-INT  PF 

—  capability 

Increase  track  lifetimes  by  developing  fusion 
techniques  resilient  to  evasive  maneuvers 


Figure  22.  TENT  Objectives 


The  Multi-INT  Particle  Filtering  Minigrant  was  performed  in-house  in  conjunction  with  Syracuse 
University.  The  final  report  is  included  as  an  appendix  to  this  report.  A  quad  chart  for  that  effort  is 
shown  in  Figure  23. 
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Multi-INT  Particle  Filtering 
Minigrant 


T+i 


Literature  Search 

Investigation  (Conceptualization) 

PF  Simulation  Development 

Analysis 

Experimentation 

Final  Report  _ 

Funding 


FY06 

/V 

- A 

A=A 

/y 

$25K 

Description:  This  project  compared  EKF  &  PF  in  the 
presence  of  glint  noise.  Changes  in  the  aspect 
toward  the  radar  can  cause  irregular 
electromagnetic  wave  reflections,  resulting  in 
significant  variation  of  radar  reflections.  This 
phenomenon  gives  rise  to  outliers  in  angle  tracking, 
and  it  is  referred  to  as  target  glint.  We  adopt  a 
commonly  used  model  for  glint  noise.  This  model 
consists  of  one  Gaussian  with  high  probability  and 
small  variance  and  another  with  small  probability 
and  very  high  variance. 

Technologies:  Multi-INT  Particle  Filtering 


Benefits  to  the  War  Fighter: 

•  Nonlinear  tracking  in  the  presence  of  glint 

•  ELINT  IMINT  fusion  using  particle  filtering 

•  Multiple  ESM  Sensor  fusion 

•  Established  In-House  Particle  Filtering  simulation 
capability 

•  Multiple  Heterogeneous  sensor  simulations 

•  Range  Bearing 

•  Bearings  Only 


•  Range  Only 


•  Established  a  baseline  approach  fortesting 
particle  filtering  algorithms  in-house 


Figure  23.  Minigrant  Quad  Chart 


The  purpose  of  the  Multi-INT  Particle  Filtering  Mini-Grant  was  to  establish  a  means  of  fusing  data  from 
Multiple  Intelligence  sources  using  Particle  Filtering  as  a  basis.  Traditional  methods  of  data  fusion 
involved  extensive  use  of  the  Kalman  Filter.  The  basic  Kalman  Filter  made  the  assumptions  that  the 
underlying  processes  were  linear  and  Gaussian.  For  example,  Electronic  Intelligence  (ELINT)  and 
Imagery  Intelligence  (IMINT)  information  require  consideration  of  the  probabilistic  characteristics  of  the 
underlying  sources  of  those  types  of  information.  It  can  safely  be  said  that  the  underlying  processes  are 
nonlinear  and  non-Gaussian.  Particle  Filters  were  typically  focused  on  the  single  information  source 
tracking  problem.  No  attempt  had  been  made  to  combine  the  state  estimation  capabilities  of  Multi-INT 
information  using  Particle  Filters.  Investigation  of  the  feasibility  of  using  Particle  Filtering  as  a  basis  for 
Multi-INT  fusion  was  accomplished  by  developing  a  MATLAB  simulation  and  experimenting  with 
promising  sensor  configurations.  Progress  was  made  in  a  number  of  areas  including: 
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1.  Established  an  in-house  software  simulation  capability  to  do  particle  filtering 

2.  Established  a  baseline  approach  for  testing  particle  filtering  algorithms  in-house 


Several  conclusions  have  been  reached  as  a  result  of  FTNLF  effort  thus  far.  The  problem  of  non-linear 
non-Gaussian  behavior  is  a  difficult  problem  that  cannot  be  addressed  through  conventional  methods. 
Investigation  of  Particle  Filtering  and  other  recent  innovations  in  non-linear  filtering  has  show  promise. 
Considerable  study  is  necessary  to  reach  conclusive  results  based  on  algorithm  performance.  Current  and 
future  work  is  focusing  on  developing  sufficient  tools  and  understanding  of  the  mathematical  nature  of 
non-linear  filtering  techniques  to  fully  determine  the  conditions  under  which  these  techniques  will  lead  to 
improved  performance. 

This  work  has  also  pointed  to  some  future  directions.  One  of  the  key  areas  that  needs  to  be  investigated  is 
the  application  of  Measures  (or  Degrees)  of  Nonlinearity  (MoN,  DoN)  to  determine  when  it  make  sense 
to  use  a  specific  type  of  nonlinear  filter  (e.g.  PF).  Also  required  is  the  investigation  into  specific  types  of 
scenarios  that  would  dictate  the  use  of  new  nonlinear  filtering  techniques.  Once  this  is  established,  the 
next  step  is  implementation  of  a  wide  variety  of  nonlinear  filtering  techniques  within  a  software  suite  for 
further  investigation  and  experimentation.  Finally,  extensive  documentation  of  results  and  conclusions 
from  using  the  experimental  suite  is  required. 
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Appendix  I.  Multi-INT  Particle  Filtering  Minigrant  report 


Particle  Filtering  Experimentation 


Summary 

This  project  develops  several  data  fusion  and  target  tracking  algorithms  for  a  surveillance  system  that 
consists  of  multiple  heterogeneous  sensors.  Algorithms  based  on  the  classical  extended  Kalman  filter 
(EKF)  and  on  the  emerging  non-Gaussian  and  nonlinear  particle  filtering  (PF)  techniques  have  been 
implemented.  These  algorithms  have  been  tested  in  the  practical  case  where  a  target  maneuvers  from  time 
to  time  and  an  Interacting  Multiple  Model  (IMM)  framework  is  used.  We  have  also  tested  them  in  the 
presence  of  glint  noise.  Changes  in  the  aspect  toward  the  radar  can  cause  irregular  electromagnetic  wave 
reflections,  resulting  in  significant  variation  of  radar  reflections.  This  phenomenon  gives  rise  to  outliers  in 
angle  tracking,  and  it  is  referred  to  as  target  glint.  We  adopt  a  commonly  used  model  for  glint  noise,  the 
Gaussian  mixture  model.  This  model  consists  of  one  Gaussian  with  high  probability  and  small  variance 
and  another  with  small  probability  of  occurrence  and  very  high  variance. 

The  performances  of  the  EKF  and  the  particle  filter  have  been  compared  through  extensive  simulation 
experiments.  The  results  show  that  for  highly  non-linear  measurements,  such  as  those  from  multiple 
bearing-only  sensors,  particle  filter  exhibits  a  superior  data  fusion  and  tracking  performance  than  the 
EKF.  However,  if  the  system  receives  measurements  from  a  radar  (both  bearing  and  range 
measurements),  the  EKF  and  the  PF  have  very  similar  tracking  accuracy,  and  the  EKF  is  a  more  desirable 
choice,  considering  that  it  requires  much  less  computation  than  the  PF,  and  has  a  much  easier  real-time 
implementation. 

Index  Terms 

Target  tracking,  Particle  Filter,  Extended  Kalman  Filter  (EKF),  Interacting  Multiple  Model  (IMM),  glint 
noise,  Gaussian  Mixture  Model  (GMM),  sensor  networks 

I.  Introduction 

Combining  the  information  from  multiple  heterogeneous  sensors  can  lead  to  more  accurate  tracking 
results  than  using  a  single  sensor.  To  fuse  these  heterogeneous  and  non-linear  measurements,  there  are 
many  tracking  algorithms,  of  which  the  most  commonly  used  is  the  classical  method  called  the  extended 
Kalman  filter  (EKF)  [1],  where  the  non-linear  measurement  model  and/or  nonlinear  motion  model  are 
linearized  via  Taylor  series  expansion,  and  the  noises  are  approximately  assumed  Gaussian.  On  the  other 
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hand,  a  Monte-Carlo  simulation  based  recursive  estimation  algorithm,  the  particle  filtering  (PF)  algorithm 
[2,  3]  has  emerged  as  a  very  promising  technique  to  solve  the  non-linear  and  non-Gaussian  filtering 
problem.  It  has  been  shown  that  using  highly  nonlinear  measurements,  such  as  bearing-only 
measurements,  the  PF  outperforms  the  EKF  [2]. 

Here  in  this  project,  we  compare  the  tracking  performance  of  the  EKF  and  the  PF  under  various 
situations,  where  different  combinations  of  sensor  measurements  are  available  for  data  fusion.  Different 
types  of  sensors  are  considered,  including  range-only  sensors,  bearing-only  sensors,  and  radars  that 
provide  both  range  and  bearing  measurements.  Besides  the  non-linearity  in  the  measurements,  non- 
Gaussian  measurement  noise,  the  glint  noise  modeled  as  a  Gaussian  mixture,  has  been  used  in  the 
experiments.  In  addition  to  the  relatively  easy  case  where  the  target  moves  at  nearly  a  constant  velocity, 
we  investigate  the  difficult  case  where  targets  maneuvers  and  an  IMM  algorithm  has  to  be  used.  Through 
simulation  experiments,  we  demonstrate  that  the  particle  filter  has  superior  performance  during  the  first 
several  steps  after  initialization.  In  steady  state,  when  the  data  are  highly  nonlinear  bearing-only 
measurements,  the  PF  still  outperforms  the  EKF.  However,  whenever  radar  data  are  available,  the  PF  has 
very  similar  steady-state  performance  as  the  EKF  in  terms  of  MSE. 

II.  Sensor  Network  Setup 

Since  multiple  heterogeneous  sensors  are  connected  to  form  a  sensor  network,  it  is  very  important  to  take 
advantage  of  the  information  from  multiple  sources.  Here  we  adopt  one  of  the  most  common  data  fusion 
schemes,  namely  the  centralized  fusion  scheme.  In  a  centralized  fusion  process,  all  the  sensors  transmit 
their  raw  measurements,  such  as  range,  and  bearing  to  the  fusion  center,  as  shown  in  Fig.  1.  After 
collecting  all  these  measurements,  the  fusion  center  fuses  them  to  form  a  new  and  more  accurate  estimate 
of  the  target  state.  The  fusion  is  accomplished  by  the  tracker  in  a  very  natural  way.  Namely  all  the  raw 
measurements  and  their  associated  accuracies  are  used  to  update  the  target  state.  The  tracker  only  needs  to 
adjust  its  measurements  equation  to  reflect  that  measurements  are  from  multiple  heterogeneous  sensors. 
The  centralized  fusion  scheme  is  optimal  in  the  sense  that  no  information  is  lost  during  the  fusion  process, 
since  the  unprocessed  raw  measurements  are  transmitted  to  the  fusion  center. 
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Fig.  1  Centralized  fusion  process. 

III.  System  Models 

A.  Target  Motion  Model 

The  maneuvering  target  motion  is  modeled  by  three  switching  dynamics  models  whose  evolution  follows 
a  Markov  chain,  also  called  a  Jump  Markov  System  (JMS)  [4,  7].  We  assume  that  at  any  time,  the  target 
moves  according  to  one  of  s  =  3  dynamic  behavior  models:  (a)  Constant  Velocity  (CV)  motion  model,  (b) 
clockwise  Coordinated  Turn  (CT)  model,  and  (c)  anticlockwise  CT  model.  Let  S=  {1,  2,  3}  denotes  the 
set  of  three  models  for  the  dynamic  motion.  Then,  the  target  dynamics  can  be  written  as 

s,-=/"‘W,)  +  vl 

where  xk  is  the  state  vector  defined  by  x^  =  [x  x  y  y] ,  k  denotes  the  discrete  time  index,  and  rke  S  is 
the  regime  variable  taking  effect  in  the  time  interval  (k-1,  k],  with  transition 

A 

probabilities^.  =Pr {rk+x  =  j\rk—  /},(/,  j  e  S )  ,  such  that  >  0,^;r..  =  1 .  The  initial  probabilities 

j 

A 

are  denoted  by  ni  =  Pr{r0  =  /}  for  i  e  S  ,  and  n .  >  0  ,  =  1 .  vk  denotes  the  white  Gaussian  noise 

i 

with  covariance  matrix  Q: 

"l/3  r3  HIT2  0  0 

HIT2  T  0  0 

Q  =  q 

0  0  1/3  J3  HIT2 

0  0  HIT 2  T 

q  is  a  scalar,  and  T  is  the  sampling  time.  For  the  CV  motion  model,  the  function  can  be  replaced 

by  the  transition  matrix  .  When  rt=l,  /•  ' ' '  (■)  corresponds  to  the  standard  CV  model 
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1  T  0  0 

0  10  0 
0  0  17 

0  0  0  1 

rk=2,3  correspond  to  clockwise  and  anticlockwise  CT  motions,  respectively. 

sin(r4;)7)  l-cos(<y^)7) 

0  cos (co[j)T)  0  -sm(co{kj)T)  2  3 

l-cos(47>0  sin  {co{kj)T) 

co[j)  co[j) 

0  sin(^7)r)  0  cos(^7)r) 

Here  the  mode-conditioned  turning  rates  are  given  by 


where  am  is  the  constant  maneuver  acceleration  parameter. 

B.  Sensor  Measurement  Model 

Three  types  of  sensors  are  used  in  our  work.  These  are  1)  ESM  sensor  that  reports  bearing-only 
measurements,  2)  range  sensor  that  reports  range  measurements  and  3)  2D  RADAR  sensor  that  reports 
range-bearing  measurements  [8].  The  measurement  model  can  be  mathematically  written  as 

z*=^(0(x*)  +  w* 

where  zJk  is  the  measurement  from  sensor  j  .  /z(0(*)  corresponds  to  three  types  of  sensor  measurement 
models,  i  =  1,2,3 

hm(xk)  =  tan1  — — 

{xk-x'  J 
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hi3\xk)  = 


tan 


f  Sj 

h  zZ 

\xk  -xJ  j 


y[iyk  ~/')2  +(xk  -x’J)2 


where  (xk ,  yk )  is  the  target  position  at  time  k,  (xSj  ,ySj)  is  the  position  of  sensor  j .  denotes  the 

measurement  noise.  In  our  work,  we  examine  the  glint  noise  as  well  as  standard  Gaussian  noise  [5].  Glint 
noise  has  a  non-Gaussian  distribution,  and  a  mixture  approach  is  widely  used  in  modeling  the  non- 
Gaussian  glint  noise.  In  the  proposed  tracking  algorithm,  the  glint  noise  is  modeled  by  a  Gaussian 
Mixture  Model  (GMM)  with  two  components. 

where  ag>  0.5  is  the  glint  probability,  and  Z!<Z2  .  Note  that  when  ag  =  1 ,  glint  noise  degenerates  to 
standard  Gaussian  noise  with  zero  mean  and  covariance  matrix  21 . 


IV.  Tracking  Algorithms 

This  section  describes  the  recursive  algorithms  implemented  for  tracking  a  single  target  using  EKF  or 
particle  filter  techniques.  Two  of  the  algorithms  are  EKF-based  and  the  other  two  are  PF-based  schemes. 
The  algorithms  considered  are  (i)  EKF-IMM,  (ii)  PF-IMM,  (iii)  EKF-Glint  Noise,  (iv)PF-Glint  Noise.  All 
four  algorithms  are  applicable  to  both  single-sensor  and  multi-sensor  scenarios. 


A.  Extended  Kalman  Filter 


Extended  Kalman  filter  is  a  minimum  mean  square  error  (MMSE)  estimator  based  on  the  Taylor  series 
expansion  [1].  The  mean  and  covariance  Pk  of  the  Gaussian  approximation  to  the  posterior 
distribution  of  the  states  can  be  derived  as  follows: 


Xk\k-\ 

0s 

£ 

II 

Pk\k-\ 

=^-,^r+a 

Kk 

+  KT 

V 

~  xk\k-\  +  Kk  (yk  h(xk\k-i  j  ^)) 

Pk 

Pk\k~  1  FkHkPk^_y 

where  Kk  is  the  Kalman  gain,  Jacobians  of  the  process  model  and  measurement  model  are  given  by 


Fk=- 


dx. 


(Vt=x*-i) 


=  FWM 
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Hk  = 


Adh(xk) 


dx, 


(xt=x*it-i) 

yk~ys 


o 


O'*-/)  +(**-**)  O'*-/)  +(**-**) 


xk-x 


0 


yk~y 


VO'*  -/)2  +(**  -V)2  yliyk  -/)2  +(**  -V)2 


For  centralized  measurement  EKF  fusion 


T  -iT 


Rk  = 


R, 


k,sx 


0  ••• 


where  Hk  s  is  the  Jacobian  of  the  measurement  model  for  each  sensor,  and  Rk  s  is  the  covariance  of  the 
measurements  model. 

For  the  maneuvering  target  tracking  problem,  the  IMM  algorithm  has  been  shown  to  be  one  of  the  most 
cost  effective  and  simple  approaches.  At  each  calculation  cycle,  the  IMM  consists  of  three  major  steps: 
interaction  (mixing),  filtering  and  combination.  At  each  time,  the  initial  condition  for  the  filter  matched  to 
a  certain  mode  is  obtained  by  mixing  the  state  estimates  of  all  the  filters  at  the  previous  time  under  the 
assumption  that  this  particular  mode  is  in  effect  at  the  current  time.  This  is  followed  by  a  regular  filtering 
step,  performed  in  parallel  for  each  mode.  Then  a  combination  of  the  updated  state  estimates  of  all  the 
filters  yields  the  state  estimate. 


B.  Particle  Filtering 

Particle  filters  represent  the  state  probability  density  function  approximately  through  a  set  of  samples  and 
implement  Bayesian  recursion  directly  on  the  samples  instead  of  dealing  with  the  exact  analytical 
functional  representations  of  the  distributions  [2,  3].  Tracking  framework  based  on  particle  filtering  will 
show  better  performance  on  nonlinear/non-Gaussian  problems.  The  recursive  Bayesian  filtering  paradigm 
provides  the  a  posteriori  PDF  p{*k  \  zn)  via  the  prediction  and  update  recursions. 

Prediction: 

P(*k  \\kP  =  | P(*k  K_i MV-l  |Zl:*-l/(V-l) 

Update: 
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p^k  K)= 


p(  zk\xMxkKJ 


p(Zj  z  ) 

where  the  state  evolution  is  described  in  terms  of  the  transition  probability: 

p(xk  I  Vi) = I  p(xk  K-,  >y*-i)/>(y*- 1  \xk-x)dyk-i 

And  how  the  given  fits  the  available  measurement  zk  is  described  as: 

P(zk  |  V )  =  {  $(zk  -  hi-  ( V » wt  ))P0 )dwk 

For  maneuvering  target  tracking,  the  aim  of  the  optimal  filter  is  to  sequentially  estimate  the  unknown 
hybrid  hidden  state  {x^r^}  given  the  observations  {zVk} .  Applying  Bayes’  rule,  the  formulation  of  the 

recursion  that  updates  ri:£-i  I  zi±-i)  to  P(xw>r\:k  I  zv.k)  can  be  derived  as 

P(Zk  \X0:k->ri:k’Zl:k-l)f(Xk  \Xk-\iVk-\)71  rk_xrk 


P(X0.kP\:k  I  Zl±)  =  P(X0±-VrX±-X  I  ZV.k-l)- 


c 


where  C  is  a  constant.  The  basic  idea  for  solving  maneuvering  target  tracking  using  a  particle  filter  is  to 
decouple  the  hybrid  estimation  problem  into  a  discrete  part  and  a  continuous  part.  We  assume  that  sensor 
measurements  are  independent  from  each  other.  Here  is  the  summary  of  the  particle  filter  solution  for  the 
maneuvering  target  tracking  problem  in  sensor  networks. 

•  Generate  samples  of  rkl)  from  an  importance  proposal  distribution  rkl)  ~  7t{rk  |  Zvk,r^k_x} , 
where  Zvk  represents  the  measurements  from  all  the  sensors  up  to  time  k.  Generate 


samples  ~  p(xk  |  xt_„/£,J)  . 
Evaluate  the  importance  weights 


X  xnf_x 


P(ff  I  ffi) 


Normalize  the  weights 


YIp(4  izi 


(0  - . 


n 7 , 


(0 


m, 


U) 


J\:k- 1 


•  Resampling:  multiply/discard  particles  {rk\i  =  1,2,..., N) with  respect  to  high/low  normalized 
importance  weights  to  obtain  N  samples  {rkl\^k}*r=l 

N 

.  MMSE  ofx,=X 

i= 1 
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V.  Simulation  Results 


We  address  the  problem  of  tracking  a  maneuvering  target  in  noise  using  multiple  heterogeneous  sensors. 
Fig  2  shows  the  tracking  scenario.  In  our  experiments,  we  use  the  dynamic  state  space  model  to  generate 
the  synthetic  data,  and  50  Monte  Carlo  computation  simulations  were  carried  out  to  evaluate  the 
performance  of  the  algorithms  for  each  experiment.  The  position  mean  square  error  is  defined  as 
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Fig  2.  Simulated  path  of  the  target 
We  set  sampling  rate  T=l,  q=20,  then 
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The  typical  maneuver  acceleration  parameter  is  set  to  am=  lm/s2.  Mode  probability  transition  matrix 
used  in  the  IMM  is 


n  = 


0.9 

0.05 

0.05 

0.6 

0.3 

0.1 

0.6 

0.1 

0.3 

We  set  the  glint  probability^  =  0.9  ,  measurement  noise  covariance  S2  =  1 021 ,  and 


2,= 


3.0462e-004  0 
0  5 


A.  Tracking  target  using  two  bearing-only  sensors 
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Bearing-only  sensors  are  located  at  (xSl  ,ySl )  =  (100m, 100m)  and  (x'2 ,  ySl )  =  (0m,250m) , 
respectively.  PF  shows  better  tracking  performance  than  EKF.  As  we  can  see  later,  this  is  the  most 
difficult  case  to  track  the  target,  and  the  tracking  results  of  using  two  bearing-only  sensors  are  much 
worse  than  those  of  using  two  range  sensors.  For  the  difficult  case  (A.2)  where  the  target  is  maneuvering, 
we  can  see  that  the  MSE  for  both  EKF  and  PF  are  higher  than  those  in  the  case  where  target  motion 
follows  a  CV  model  (A.l). 


A.l  CV  model  with  glint  measurement  noise 
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B.  Tracking  target  using  one  radar  sensor. 

We  assume  that  a  radar  is  located  at  ( x\ys )  =  (500m,  500m) .  From  the  experimental  results, 
we  can  see  that  except  for  the  first  few  steps,  EKF  and  PF  achieve  almost  the  same  MSE.  But  the 
computation  time  is  much  shorter  for  the  EKF.  For  the  difficult  case  where  the  target  is  maneuvering,  we 
can  see  that  the  MSE  for  both  EKF  and  PF  are  higher  than  those  for  the  case  of  target  motion  that  follows 
a  CV  model. 

B.l  CV  model  with  glint  measurement  noise 
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EKF  Estimate  &  PF  vs  True 
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B.2  Maneuvering  target 
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C.  Tracking  target  using  one  range  sensor  and  one  bearing-only  sensor 

The  range  sensor  position  is  set  to  (xSl  ,ySl )  =  ( 0m,250m )  ,  and  bearing-only  sensor  position  is 
set  to  (xSl  ,ySl )  =  (100m, 100m)  .  The  function  of  bearing-only  sensor  plus  range  sensor  is  almost  the 
same  as  a  single  radar  sensor,  except  the  bearing-only  sensor  and  range  sensor  are  located  at  different 
locations,  so  we  have  similar  experimental  results  as  in  case  B. 

C.l  CV  model  with  glint  measurement  noise 
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C.2  Maneuvering  target 
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EKF  Estimate  &  PF  vs  True  Ekf  &  PF  MSE 


D.  Tracking  target  using  two  range  sensors 

Two  range  sensors’  positions  are  set  to  (xs‘ , ys' )  =  (100m, 250m)  and  (x*2 , )  =  (300m, Om) , 
respectively.  For  this  sensor  configuration,  the  PF  still  has  similar  steady-state  performance  as  that  of 
the  EKF. 

D.l  CV  model  with  glint  measurement  noise 
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E.  Tracking  target  using  bearing  only,  range  and  radar  sensors 

Sensors’  positions  are  set  to(x*' ,  jSl )  =  (250m, 0m),  (x*2 ,  j*2 )  =  (0m,250m) 

and(x*3  ,yS} )  =  (500m, 500m) .  Here  we  use  three  different  sensors,  and  MSE  is  much  smaller  than  the 
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previous  cases.  As  expected,  more  accurate  tracking  results  are  achieved,  since  we  are  fusing  data  from 
more  sources.  Again,  the  PF  and  the  EKF  have  very  close  steady-state  performance. 

E.l  CV  model  with  glint  measurement  noise 


EKF  Estimate  &  PF  vs  True 
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From  the  above  experiments,  we  can  observe  that  even  when  the  measurement  model  is  nonlinear  for  all 
types  of  sensors,  using  particle  filtering  does  not  always  achieve  a  better  steady-state  performance  in 
terms  of  MSE.  Only  under  the  conditions  where  only  bearing-only  sensors  are  used,  the  particle  filter 
shows  better  performance  than  the  EKF.  We  also  found  that  even  when  the  initial  condition  for  both  EKF 
and  PF  is  the  same,  in  the  first  few  tracking  steps,  PF  tracking  results  are  more  accurate  than  EKF.  This  is 
a  valuable  characteristic,  especially  when  clutter  and  false  alarms  are  among  the  measurements.  When  the 
measurements  contain  many  false  alarms,  there  is  uncertainty  as  to  which  measurement  is  from  the  target 
and  which  is  a  false  alarm.  With  such  uncertainty,  inaccurate  estimates  even  at  one  time  step  could  lead 
the  filter  to  diverge  and  result  in  the  loss  of  the  target  track.  The  PF  has  the  potential  to  maintain  the  target 
track  for  a  longer  time  in  such  harsh  and  realistic  conditions.  This  issue  needs  further  investigation  in  the 
future. 
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Introduction 


During  the  summer  of  2007  I  worked  in  the  Fusion  Technology  branch  (IFEA).  I 
was  part  of  an  in-house  team  focused  on  developing  and  improving  non-linear  tracking 
techniques.  I  aided  in  the  finalization  to  the  Development  and  Evaluation  of  Fusion 
Techniques  program  known  as  DEFT.  I  also  worked  on  the  In-House  Particle  Filtering 
and  Testing  program  along  with  the  newest  effort:  Tracking  Evasive  Non-Linear  Targets 
(TENT). 

These  efforts  in  one  way  or  another  were  related  to  non-linear  tracking.  I  focused 
on  developing  and  evaluating  different  techniques  to  deal  with  non-linear  tracking.  I 
researched  a  variety  of  filters  based  upon  the  Kalman  filter  and  the  particle  filter.  Other 
topics  of  interest  I’ve  examined  include:  out  of  sequence  measurements  (OOSMs), 
multiple  target  tracking,  flight  dynamics,  sensor  characteristics,  methods  for  evaluating 
tracker  performance,  probability  theory  and  statistics,  as  well  as  many  others. 

In  this  report  a  basic  introduction  to  particle  filters  illustrating  the  details  of  its 
operation  as  well  as  its  pros  and  cons  are  provided.  Thorough  analysis  of  in-house 
simulations  run  during  the  summer  is  presented.  Two  sources  for  code  are  presented  for 
future  reference.  Finally  recommendations  for  future  research  are  detailed. 

Background  on  Particle  Filters 

Kalman  filters  have  been  used  for  decades  and  various  adaptations  have  been 
produced  through  the  years  to  deal  with  troubling  issues  such  as  dealing  with  non¬ 
linearity.  A  solution  used  to  address  some  of  the  weaknesses  exhibited  in  Kalman  filters 
is  the  particle  filter.  Although  particle  filtering  has  been  around  for  quite  some  time  too, 
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it  has  only  recently  come  back  into  the  spotlight.  In  the  past  the  particle  filtering  was 
deemed  too  computationally  expensive  for  realistic  use  but  with  the  advent  of  faster 
processors  and  technological  upgrades,  particle  fdters  are  acquiring  a  new  look. 

“Particle  fdters  are  sequential  Monte  Carlo  methods  based  upon  point  mass  (or 
‘particle’)  representations  of  probability  densities,  which  can  be  applied  to  any  state 
space  model,  and  which  generalize  the  traditional  Kalman  filtering  methods”  [1],  In 
other  words,  a  particle  fdter  can  estimate  a  state  space  model  by  representing  the 
posterior  density  function  using  many  random  samples  (particles)  with  assigned  weights. 
As  the  number  of  these  particles  increases  the  representation  of  the  pdf  reaches  its 
optimal  state. 

There  are  numerous  variations  of  particle  fdters  but  one  of  the  most  basic  is  a 
Sequential  Importance  Resampling  (SIR)  Filter.  There  are  only  a  few  basic  steps  to  this 
fdter.  Initialization  begins  with  establishing  a  prior  probability  distribution,  that  is,  an 
initial  guess  in  the  form  of  a  probability  distribution  over  the  state  at  the  start  time.  This 
corresponds  to  how  likely  it  is  for  the  target  to  be  at  a  given  location.  Next  a  set  of 
particles  are  drawn  from  this  initial  prior  probability  distribution  and  the  filter  is 
initialized. 

Once  the  fdter  is  initialized  it  repeats  three  basic  steps.  The  first  is  generating  a 
proposal  distribution.  This  step  asks  “where  could  the  target  have  moved  to  given  what  I 
knew  about  where  it  might  have  been  an  instant  before.”  It  predicts  the  state  of  the  target 
given  all  previous  observations. 
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The  next  step  is  incorporating  observations.  This  step  takes  information  from 
the  sensors  and  updates  the  belief  where  the  target  is  located  by  assigning  importance 
weights  to  all  the  particles. 

The  final  step  is  resampling  and  it  is  required  for  the  filter  to  maintain  diversity. 
SIR  filtering  is  unique  in  that  the  resampling  step  aids  in  reducing  the  degeneracy 
problem.  After  a  few  iterations  of  the  basic  particle  filtering  process  all  but  one  particle 
will  have  negligible  weight.  Resampling  determines  the  number  of  effective  particles 
(i.e.  particles  with  high  weights)  remaining  and  if  it  is  below  a  certain  threshold  shifts  the 
negligible  particles  closer  to  the  effective  particles.  The  filter  resamples  particles 
according  to  the  renormalized  weights.  Particles  with  higher  weights  get  sampled  more 
often  than  particles  with  lower  weights. 

The  particle  filter  has  several  benefits  compared  to  other  filters.  The  particle  filter 
performs  well  in  non-linear  situations  and  is  capable  of  handling  a  non-linear  state  and 
observation  model.  The  particle  filter  is  also  designed  to  perform  well  with  non-Gaussian 
distributions  and  allows  the  use  of  multi-modal  distributions.  Unlike  some  filters  the 
particle  filter  estimates  the  full  probability  density  function  of  the  state.  With  a  near 
infinite  number  of  particles  the  results  produced  are  close  to  the  optimal  solutions. 

There  are  also  several  drawbacks  with  using  a  particle  filter.  The  most  common 
concern  is  its  computational  cost.  As  the  particles  increase  in  number,  the  number  of 
computations  increases  dramatically.  Along  with  other  filters  the  particle  filter  is  also 
prone  to  the  curse  of  dimensionality.  As  the  number  of  dimensions  increase  the 
computational  complexity  increases  exponentially.  The  problems  associated  with 
degeneracy  mentioned  previously  exist  in  many  versions  of  particle  filters.  Often  particle 
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filters  simply  aren’t  necessary  in  some  applications  and  a  basic  linear-Gaussian 
assumption  is  sufficient  to  produce  the  desired  results. 

In-House  Simulations 

Utilizing  in-house  software  already  developed  in  collaboration  with  Syracuse 
University,  I  compared  the  effectiveness  of  an  Extended  Kalman  filter  (EKF)  against  a 
particle  filter  (PF).  The  simulations  run  via  MATLAB  tracked  a  moving  ground  target 
using  an  IMM-EKF  and  a  MM-PF  with  bearings-only  sensors.  The  path  of  the  target 
could  be  modeled  in  a  variety  of  ways  to  increase  or  decrease  the  degrees  of  linearity. 
Any  number  of  bearings-only  sensors  could  be  placed  around  the  target.  I  analyzed  the 
simulation  and  determined  how  numerous  parameters  interacted  with  one  another. 

Comparison  of  the  EKF  vs.  PF  depends  on  several  factors  including:  sensor 
position,  process  noise  (determined  by  a  scalar  multiplier  “q”),  number  of  particles  used 
in  the  particle  filter,  measurement  noise,  and  glint  noise. 

The  first  scenario  that  was  extensively  analyzed  consisted  of  the  ground  target 
initially  heading  160  degrees  from  North  at  70  mph.  After  20  seconds,  the  target 
executed  a  Counterclockwise  Coordinated  Turn  for  3  seconds,  to  establish  a  new  course. 
The  target  then  traveled  with  a  Constant  Velocity  motion  until  it  executed  its  final 
maneuver  for  4  seconds,  a  Clockwise  Coordinated  Turn. 

For  this  scenario,  two  bearings-only  sensors  were  activated  at  a  time  and 
positioned  at  fixed  points  surrounding  the  moving  target.  To  understand  the  effects  of 
sensor  positioning,  these  two  sensors  were  repositioned  in  15  different  setups. 
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Sensor  position  can  degrade  overall  tracking  performance  for  both  EKF  and  PF. 
Overall  tracking  performance  refers  to  RMSE  and  variance.  A  high  tracking 
performance  is  described  as  having  low  RMSE  and  low  variance.  Performance  is 
degraded  when  the  sensors  are  not  lined  up  normal  to  the  target  track.  When  the  sensors 
are  lined  up  in  parallel  with  the  direction  of  motion  of  the  target,  performance  is 
minimized. 

In  this  scenario,  overall  tracking  performance  was  at  a  minimum  when:  1)  the 
two  sensors  were  positioned  side  by  side  on  all  four  sides  of  the  target’s  track;  2)  the  two 
sensors  were  positioned  on  opposite  sides  of  the  target  track  parallel  with  the  target’s 
motion.  Alternatively,  tracking  performance  was  maximized  when  either:  1)  one  or 
more  of  the  sensors  were  positioned  close  to  the  target  track  so  at  no  time  during  the 
simulation  was  the  sensor  in  parallel  with  the  target’s  motion;  2)  the  two  sensors  were 
positioned  at  right  angles  on  the  outside  of  the  target  tack  ensuring  at  least  one  sensor  was 
normal  to  the  target  at  all  times. 

In  the  two  cases  mentioned  previously  where  performance  is  minimized,  the  EKF 
outperformed  the  PF  in  nearly  every  situation.  When  the  bearings-only  sensors  are 
positioned  in  such  a  manner  the  particle  fdter  has  difficulty  recognizing  range  and  as  a 
result  may  make  a  wrong  turn,  while  the  EKF  does  not  experience  this  problem.  This 
may  be  due  to  the  specific  implementation  of  each  filter  or  perhaps  may  suggest  the 
particle  filter  is  more  suitable  with  an  addition  of  a  range  sensor  such  as  radar.  Despite 
the  EKF’s  performance  over  the  PF,  it  is  important  to  remember  the  RMSE  and  the 
variance  were  still  very  large  for  both  filters  in  this  sensor  setup. 
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When  the  sensors  are  positioned  in  the  two  cases  where  performance  is 
maximized  the  results  weren’t  as  one-sided  toward  the  EKF  as  when  the  sensors  were 
positioned  poorly.  A  noticeable  trend  emerged  from  this  sensor  setup.  The  PF 
outperformed  the  EKF  whenever  “q”  was  set  to  a  negligibly  low  value  no  matter  how 
many  particles  were  used  (50,  100,  500).  This  poor  performance  by  the  EKF  with  a  low 
“q”  value  was  found  to  be  the  direct  result  of  an  initialization  issue.  In  this  scenario  the 
two  filters  initialized  their  tracks  at  a  different  position  than  the  actual  initial  target 
location.  This  caused  the  EKF  to  lose  the  track  immediately  after  the  first  few  time  steps 
because  the  low  “q”  value  represented  high  confidence  in  the  EKF’s  measurements  when 
in  fact  the  measurement  was  far  away  from  the  actual  target  location.  The  particle  filter 
did  not  exhibit  any  such  behavior  with  the  varying  process  noise  levels.  Thus  it  can  be 
seen  that  the  EKF  is  highly  dependent  on  process  noise  while  the  particle  filter  is  less 
dependent. 

When  the  PF  was  set  to  50  particles,  in  this  performance  maximization  sensor 
setup  once  again,  the  EKF  showed  better  performance  for  the  “q”  values  of  0.5,  1,  and  5. 
As  the  number  of  particles  increased,  the  performance  of  the  PF  increased  (with  the 
exception  that  simulation  time  increased  dramatically  as  a  result  from  more 
computational  complexity).  When  the  number  of  particles  was  increased  to  100  the  EKF 
outperformed  or  did  just  as  well  as  the  PF  with  the  exception  of  when  “q”  was  set  to  a 
negligibly  small  number  such  as  1.6  xl0~6  as  mentioned  above.  Once  the  number  of 
particles  was  set  to  500  the  PF  either  outperformed  or  did  just  as  well  as  the  EKF.  As 
expected  theoretically,  as  the  number  of  particles  approaches  infinity,  the  closer  the 
results  will  become  optimal. 
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The  effects  of  varying  levels  of  measurement  noise  were  briefly  evaluated.  For 
ranging  levels  of  measurement  noise  it  was  discovered  there  exists  a  maximum  and 
minimum  threshold  not  to  be  exceeded.  Performance  deteriorates  severely  when  the 
measurement  noise  becomes  too  high  or  even  too  low. 

The  ensuing  series  of  scenarios  generated  focused  on  the  degrees  of  non-linearity 
of  the  target’s  track.  The  main  independent  variable  was  track  configuration.  Four 
different  paths  for  the  ground  target  were  examined:  a  straight  line,  a  circle,  a  u-tum,  and 
a  multiple  turn  model.  The  parameters  tested  in  the  previous  scenario:  process  noise, 
measurement  noise,  and  number  of  particles  were  set  to  fixed  values.  Multiple  bearings- 
only  sensors  were  positioned  in  fixed  locations  around  all  four  target  tracks  in  an  attempt 
to  minimize  the  impending  effects  of  sensor  positioning  discovered  in  the  first  scenario. 

Results  from  these  four  different  target  paths  failed  in  their  attempt  at  providing 
insight  on  the  effects  of  target  non-linearity.  For  all  four  target  paths  the  EKF  and  PF 
produced  nearly  identical  results  and  no  visible  pattern  could  be  established  as  to  when 
one  filter  outperformed  the  other.  However  this  may  be  due  to  a  fundamental  flaw  in  the 
simulation  itself.  Sampling  from  the  measurements  occurs  at  every  time  step  and  as  a 
result  may  be  sampling  too  quickly.  With  such  a  high  measurement  sample  rate  the  two 
filters  may  be  impersonating  a  non-linear  trajectory  such  as  a  circle  as  linear. 

Although  increasing  the  number  of  sensors  in  these  scenarios  was  intended  to  be  a 
panacea  for  dealing  with  sensor  alignment  issues  discussed  in  the  first  scenario,  it  failed 
in  simplifying  the  experiments.  With  an  increase  in  the  number  of  sensors  tracking 
performance  is  naturally  increased.  However  with  increased  performance  for  both  filters 
it  becomes  difficult  to  visually  substantiate  one  filter’s  performance  over  the  other.  If  the 
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number  of  sensors  is  reduced,  performance  is  reduced  to  the  point  where  it  is  difficult  to 
determine  whether  a  specific  fdter  reacting  to  non-linearity  is  to  blame  or  whether  it’s 
just  poor  sensor  positioning  coupled  with  a  high  sampling  rate. 

The  effect  of  glint  noise  was  also  examined  for  the  multiple  turn  model  of  the 
target  trajectory.  It  was  clear  that  glint  noise  had  a  detrimental  effect  on  performance. 
As  the  effects  of  glint  noise  were  increased,  performance  decreased.  Although  the 
performance  did  not  change  immensely  as  the  effects  of  glint  were  increased  (perhaps 
due  to  the  use  of  multiple  sensors  and  high  sample  rate),  the  transient  response  of  the 
generated  tracks  showed  a  moderate  change  in  performance.  More  importantly,  and 
perhaps  of  high  note,  the  PF  initially  outperforms  the  EKF  in  the  presence  of  glint.  This 
may  prove  the  theory  that  particle  fdters  can  outperform  an  EKF  in  non-Gaussian  noise 
environments. 

Other  Examples  of  Code 

In  an  undertaking  to  fully  understand  the  efficacy  of  a  particle  filter,  past 
implementations  of  particle  filters  already  constructed  were  investigated.  Two  main 
sources  are  to  be  highlighted.  The  first  is  code  written  by  Michael  A.  Goodrich  to 
support  Brad  Huber’s  thesis  work  [2]  and  the  second  source  is  a  series  of  examples  from 
Dan  Simon’s  book  Optimal  State  Estimation:  Kalman,  H-infinity,  and  Nonlinear 
Approaches  [3], 

Huber’s  thesis  proposed  to  equip  a  mini-UAV  with  radio  sensors  and  use 
Bayesian  filtering  techniques  to  locate  an  Emergency  Locator  Transmitter  for  a  downed 
aircraft.  Huber  created  a  MATLAB  simulation  of  a  mini-UAV  with  noisy  non-linear 
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sensor  behavior  that  searches  for  a  downed  aircraft.  Goodrich  adapted  Huber’s  code  and 
created  a  simplified  version  of  the  existing  particle  filter  implementation. 

The  code  is  important  because  it  acts  as  an  analysis  tool.  One  can  directly  see  the 
effects  of  changing  particle  filter  related  parameters  on  the  simulation.  Since  most  of  the 
information  available  on  particle  filtering  is  theory  and  pseudo  code,  the  code  provides  a 
rare  opportunity  for  seeing  actual  values. 

The  code  from  Dan  Simon  provides  several  examples  of  not  only  particle  filtering 
but  other  examples  of  non-linear,  Kalman,  and  H-infinity  filtering.  Simon’s  examples  are 
basic  and  designed  to  be  simple  supplemental  guides  to  his  book.  Studying  these 
examples  provide  all  the  same  advantages  as  Huber’s  thesis  work  but  also  offer  variety 
and  alternative  techniques. 

Recommendations  for  Future  Research 

While  understanding  that  some  initial  research  may  have  already  taken  place  and 
that  it  takes  time  to  produce  results,  recommendations  for  future  research  are  still 
presented. 

When  trying  to  produce  a  simulation  it  is  absolutely  essential  that  the  confines 
and  parameters  are  known.  For  the  future,  work  needs  to  be  done  researching  “real 
world”  values  to  apply  to  simulations.  Data  needs  to  be  collected  using  existing 
implementations  of  non-linear  tracking  as  well  as  sensor  information.  Research  should 
not  only  include  reading  for  information  but  cataloging  the  results  of  what  has  been 
learned. 
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Development  of  a  sophisticated  scenario  generator  needs  to  be  developed.  A 
generator  that  is  capable  of  handling  multiple  tracking  filters  with  an  easy  to  use  interface 
for  changing  parameters  would  be  the  ultimate  goal.  The  tracks  and  values  generated 
should  be  based  on  realistic  recorded  data.  Methods  to  evaluate  tracking  performance 
such  as  root  mean  square  error,  variance,  or  Cramer-Rao  Lower  Bounds  should  be 
included. 

If  future  work  on  the  existing  software  is  desired,  changes  should  be  immediately 
made  to  the  measurement  sample  rate.  If  new  results  are  generated  from  that  change  one 
may  be  able  to  see  the  effects  from  varying  parameters  more  clearly.  Introduction  of 
alternative  types  of  sensors  such  as  radar  may  show  different  results.  Further  research 
into  the  effects  of  non-Gaussian  noise  on  different  filters  should  be  looked  into. 

It  is  necessary  to  establish  a  general  base  of  knowledge  before  advanced  measures 
can  be  introduced.  While  remembering  to  take  a  “crawl  before  you  walk”  approach, 
future  work  in  implementing  advanced  measures  should  be  implemented  after 
preliminary  results  have  been  produced.  Meaning  it  is  important  to  understand  how  the 
elements  of  non-linear  tracking  interact  with  one  another  at  the  fundamental  level  first 
and  foremost.  Once  this  is  achieved  more  complex  incorporations  may  be  introduced  to 
simulations  such  as  Out  of  Sequence  Measurements  (OOSMs),  multiple  target  tracking, 
different  filter  types,  adaptive  filters,  smoothing  techniques,  or  any  others. 
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Appendix  III.  Summer  Student  Report  2007:  Becky  Bailey 

Nonlinear  Filtering  Using  Solutions  to  the  Fokker-Planck 
Differential  Equation 

Target  tracking  has  been  one  of  the  problems  encountered  by  engineers  for  many 
decades.  The  main  goal  of  this  problem  is  to  track  a  target  that  is  moving  on  some 
trajectory  using  measurements  given  by  sensors.  Although  on  the  surface,  this  may  seem 
like  a  fairly  simple  problem,  there  are  numerous  challenges  contained  within  the  different 
tracking  scenarios.  These  include  the  facts  that  sensors  produce  noisy  measurements, 
data  is  often  received  at  irregular  time  intervals,  obstructions  occurring  in  the  “view”  of 
the  target  (i.e.  bridges,  buildings,  etc.),  data  association  problems  due  to  multiple  targets, 
and  nonlinearity  of  tracks.  Due  to  the  previously  mentioned  challenges  that  occur  in 
target  tracking,  developing  an  algorithm  that  solves  the  entire  tracking  problem  and  takes 
into  account  the  many  different  scenarios  that  can  develop  has  been  a  very  difficult  task 
for  engineers.  This  paper  will  focus  on  one  of  the  main  problems  in  tracking,  which  is 
the  nonlinearity  of  a  target’s  track.  More  specifically,  a  new  nonlinear  tracking  filter  will 
be  explained  which  uses  solutions  of  the  Fokker-Planck  differential  equation  to  solve  the 
nonlinear  tracking  problem. 

In  general,  tracking  algorithms  can  be  divided  into  three  main  steps:  initialization, 
prediction,  and  update.  All  tracking  filters  that  have  been  developed  deal  with  these  three 
steps  in  some  form  or  another.  The  initialization  step  refers  to  the  process  of  obtaining  an 
initial  state  based  on  a  certain  degree  of  uncertainty.  The  prediction  step  predicts  where 
the  target  should  be  based  only  on  previous  measurements  using  the  system  dynamical 
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model.  This  creates  the  state  estimate.  The  update  step  actually  uses  the  measurement 
upon  arrival  from  the  sensor  to  update  the  state  estimate.  A  probabilistic  based  weighting 
is  used  to  create  an  optimal  combination  between  old  and  new  information.  [1] 

The  majority  of  filtering  algorithms  developed  up  to  this  point  have  either  dealt 
with  linear  tracks,  or  used  some  sort  of  linearization  process  to  linearize  the  nonlinear 
target  tracks.  One  of  the  most  widely  known  filtering  algorithms  developed  for  target 
tracking  is  the  Kalman  Filter,  which  was  developed  in  the  1960’s  by  Rudolf  Kalman. 
Although  this  filter  is  the  optimal  filter  for  the  linear-Gaussian  tracking  problem,  its  real 
world  application  is  minimal  due  to  the  numerous  assumptions  made  in  the  theory.  In 
particular,  the  Kalman  Filter  not  only  assumes  the  target  track  to  be  linear,  but  the  noise 
signal  and  the  posterior  density  are  both  assumed  to  be  Gaussian.  These  are  very  limiting 
assumptions  because  rarely  does  this  type  of  scenario  occur  in  the  real  world,  making  it 
obvious  that  the  Kalman  Filter  needed  to  be  extended  to  nonlinear  situations.  This  was 
done  with  the  Extended  Kalman  Filter  (EKF),  which  basically  linearized  the  nonlinear 
situations  by  assuming  local  linearity  of  the  target’s  track. 

The  filters  discussed  thus  far  have  not  really  “solved”  the  nonlinear  tracking 
problem,  but  have  linearized  the  nonlinear  problem,  and  then  calculated  a  linear  solution. 
While  there  have  been  some  other  nonlinear  tracking  filters  that  do  not  actually  solve  the 
Fokker-Planck  Equation  (FPE),  such  as  the  particle  filter,  my  focus  during  my  time  here 
has  been  on  nonlinear  filters  that  solve  the  Fokker-Planck  Equation.  To  understand  these 
filters,  and  their  approach  to  the  nonlinear  tracking  problem,  one  must  first  have  an 
understanding  of  what  the  FPE  is  and  how  it  is  used  in  tracking. 
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The  Fokker-Planck  Equation  is  a  second  order  partial  differential  equation  of 
motion  for  the  probability  density  function  of  continuously  changing  macroscopic 
variables.  Solving  this  equation  leads  to  a  probability  density  function  enabling  one  to 
find  any  of  the  averages  for  macroscopic  variables  by  integration.  When  applied  to 
nonlinear  filtering,  this  equation  explains  the  progression  of  the  conditional  probability 
density  of  the  state  vector  between  the  times  at  which  measurements  are  actually  taken. 
To  solve  this  equation  in  real  time  would  solve  the  nonlinear  filtering  problem.  The 
difficulty  in  this  is  that  the  solution  suffers  from  the  “curse  of  dimensionality,”  and 
therefore  the  computational  complexity  of  the  solution  grows  exponentially  with  the 
dimension  of  the  state  vector.  [2] 

The  general  FPE  is  given  in  the  following  form: 


dW_ 

dt 


W(x,t) 

i= 1  VXi  i,j=\  OX ; 


(l) 


With  the  variables  being  defined  as  follows: 
D\V)  =  Drift  vector 

D-P  =  Diffusion  tensor 
W({x},t)  =  Distribution  function 


The  FPE  is  actually  defined  by  the  first  2  terms  of  the  Kramers-Moyal  expansion.  The 
reason  that  it  is  only  the  first  two  terms  of  the  Kramers-Moyal  expansion  is  that  the 
expansion  with  an  infinite  number  of  terms  stops  after  the  second  term  due  to  the  fact  that 
the  coefficients  vanish  for  n>3,  according  to  the  Langevin  equations  from  which  they 
are  derived.  The  drift  vector  and  diffusion  tensor  are  therefore  just  the  first  two 
coefficients  of  the  Kramers-Moyal  expansion. 
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In  nonlinear  filtering,  the  FPE  takes  a  specific  form  based  on  the  model  used  to 
define  the  dynamic  system.  This  model  is  given  by  the  following  equation: 
x'  =  f(x',t)  +  co  (2) 

Where  x'  is  the  system  state  vector,  while  the  forcing  function,  co  is  a  random  process. 
The  probability  density  function  p(x,t)  for  the  state  variables  defined  by  (2)  satisfies  the 


FPE  of  the  following  form: 


dp 

dt 


-^\f-p\tr\—  +~tr  Q 


(3) 


dp  8f 

In  this  case,  — —  is  the  gradient  of  p  with  respect  to  x' ,  —  is  the  Jacobian  of  /  with 

dx  dx 


d2  p 

respect  tox' ,  — —  is  the  Jacobian  of  the  transpose  of  p  with  respect  to  x' ,  and  Q  is  the 


process  noise  matrix.  [3] 

In  his  paper  [4],  Fred  Daum  develops  a  filter  which  solves  the  FPE,  and  in 
combination  with  an  exponential  function,  he  solves  the  nonlinear  filtering  problem. 
Unlike  the  Kalman  filter  which  can  only  handle  the  class  of  Gaussian  probability 
densities,  Daum’s  filter  generalizes  this  so  that  any  multivariate  probability  density  from 
the  exponential  family  can  be  used.  In  this  theory,  Daum  considers  the  random  variable 
x(t),  which  evolves  continuously  in  time  according  to  the  following  equation: 
dx(t)  =f{x,  t)dt  +  G(t)dw  (4) 

In  this  case, /is  a  nonlinear  function  of  x,  and  w(t)  is  an  Revalued  Brownian  motion.  He 
then  assumes  that  there  exists  a  positive  real- valued  function,  ¥=  W(x,t)  that  satisfies  the 
following  FPE  corresponding  to  (4): 
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dp 

~dt 


r  dxp^ 

f-w 

tr 

fs/Y 

1 

+  —tr 

f  rsl  \ 

QeV 

dx  J 

ISxJ 

2 

r  J 

With  the  initial  condition: 


ip(x,t0)  = 


P(x,t0)e  2 


{(h(x-m0)TP0-l(x-m0)} 


(5) 


(6) 


If  all  of  Daum’s  assumptions  and  conditions  are  satisfied,  the  result  of  his  theory  is  an 
exact  formula  for  the  unnormalized  probability  density  of  the  state  at  a  specific  time 
conditioned  on  the  set  of  discrete  time  measurements  Zk={zi,  Z2,...,  Zk}.  This 
unnormalized  probability  density  is  given  by  the  following  formula: 


(Xk  >  h  ) ' 


exp 


~mkY  pkl(xk 


(?) 


Where  mk  and  Pu  can  be  computed  recursively  using  the  following  formulas: 

mk  =mk+PkHTkR-k\zk - Hkmk ) 

n  =n-  prf  ( HtptHi + Rk  Hkpk  <8) 

Although  the  recursive  formulas  for  nit  and  Pk  lead  the  reader  to  believe  that  they 
are  the  mean  and  covariance  of  the  state,  they  are  not.  In  Remark  3  of  his  paper,  Daum 
states,  “Note  that  mk  is  not  the  conditional  mean  of  jc*,  and  Pk  is  not  the  conditional 
covariance  matrix  of  jc*.  But  rather,  nik  and  Pk  can  be  thought  of  as  the  mean  and 
covariance  matrix  of  some  Revalued  auxiliary  Gaussian  random  variable.”  [4]  Since  it  is 
very  difficult  to  create  a  recursive  filter  without  the  mean  and  covariance,  in  [3]  Schmidt 
develops  a  way  to  relate  Daum’s  variables  mk  and  Pk  to  the  actual  mean  and  covariance  x 
and  M.  The  following  equations  were  derived  by  Schmidt  for  the  mean  and  covariance 
based  on  Daum’s  theory: 


M  =  BM  +  MBt  - M(2L  +  gxGx  + ...  +  gnGn)M  +  Q 
AND 

5c  =  h  +  Bx-M(2L  +  gl  a +...  +  gnGn)x-MS 


(9) 


As  a  reader  familiar  with  tracking  can  quickly  see,  the  equation  for  propagating  the 

covariance  only  differs  from  that  of  the  EKF  by  the  factor  of 

M(2L  +  gfix  + ...  +  gnGn)M  (10) 

It  is  also  evident  that  the  equation  for  propagating  the  mean  only  differs  from  that  of  the 
EKF  by  the  factor  of 

M(2L  +  glGl+...  +  gnGn)x  (11) 

These  two  factors,  (10)  and  (11),  are  the  part  of  the  propagation  equations  that  take  into 
account  Daum’s  parameters  m*  and  Pk-  The  update  equations  are  also  very  similar  to 
those  of  the  EKF  and  are  given  as  follows: 

Mk(+y'xk(+)  =  HTRlzk  +Mk(-ylxk(-) 

AND  (12) 

Mk(+yl=HTR-lH  +  Mk(-y 

Where  the  parameters  before  update  are  indicated  by  (-),  and  after  update  they  are 
indicated  by  (+). 

Overall,  Daum  was  able  to  divide  the  tracking  estimation  problem  into  two 
separate  steps.  One  step  consists  of  the  recursive  calculation  of  nik  and  Pk  based  on  the 
sensor  measurements.  As  previously  revealed,  this  can  be  done  with  equations  very 
similar  to  those  of  the  Kalman  filter.  Schmidt  took  this  step  one  level  further,  and 
developed  propagation  and  update  equations  for  the  mean  and  covariance  based  on 
Daum’s  parameters  mu  and  Pk-  The  other  general  step  in  Daum’s  theory  consists  of  the 
calculation  of  W{x,i),  which  does  not  depend  on  the  sensor  measurements.  This 
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calculation  can  be  done  off-line,  before  the  actual  observation  interval  using  numerical 
solutions  to  the  FPE. 

Daum’s  filter  that  was  previously  described  is  a  “finite-dimensional”  filter,  which 
means  that  it  can  be  implemented  by  integrating  a  given  number  of  ordinary  differential 
equations.  [4]  In  [5],  Daum,  along  with  Lambert  and  Weatherwax,  develops  a  different 
filter,  which  they  call  the  “wave  filter.”  This  algorithm  solves  the  nonlinear  tracking 
problem  by  simply  solving  the  FPE  for  the  conditional  probability  density  function.  This 
is  done  in  a  “split-step  solution,”  which  suffers  from  the  “curse  of  dimensionality,”  unlike 
Daum’s  previously  explained  filter.  The  idea  behind  this  filter  is  to  solve  the  FPE  using 
the  “split-step  solution,”  which  can  be  viewed  as  the  propagation  step  of  the  algorithm. 
Then  for  the  update  step,  the  conditional  density  is  updated  upon  arrival  of  new 
measurements  using  Bayes’  rule.  This  is  given  by: 


(13) 


In  constructing  the  “split-step  solution,”  the  following  system  dynamic  model  was 
used  to  describe  the  dynamic  behavior  of  the  system: 

dxt  =  f(xt)  +  dwt  (14) 


In  this  model,  the  process  noise  vector  wt  is  a  zero  mean  Gaussian  white  noise  process 
with  spectral  density  matrix  Qt  [5].  The  FPE  used  to  describe  the  evolution  of  the 
conditional  probability  density  between  measurements  is  given  by: 

(15) 


dp 

dt 


"T  dp 


OX. 


rdft A 


\dxtj 


1 


'  'p  +—tr\ 
2 


a 


d2p\ 

dx?) 
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The  solution  described  here  will  consider  two  separate  special  cases  of  this  FPE  in  order 

to  come  to  one  unified  solution.  The  first  case  considered  will  be  the  case  where  there  is 

time  invariant  process  noise,  and  therefore,  f  =  0  and  Q,  =  Q.  With  these  conditions  in 

mind,  the  FPE  reduces  to:  —  =  —tr  O — 

8t  2  {  dx] )  (16) 

This  reduced  form  of  the  FPE  can  actually  be  viewed  as  the  diffusion  equation  with  Q 
being  the  diffusion  coefficient.  Equation  (15)  can  be  solved  by  computing  the 
convolution  integral  using  an  algorithm  such  as  the  Fast  Fourier  Transform  (FFT).  The 
second  special  case  considered  in  this  split-step  solution  is  the  case  where  there  is  no 
process  noise,  and  therefore  Qt  =  0.  In  this  case,  the  FPE  reduces  to: 


dp+  Tdp 
dt  1  dx, 


=  —tr\ 


P 


(17) 


In  this  form,  f  can  be  viewed  as  the  drift  vector,  while  the  equation  itself  has  actually 
been  reduced  to  the  convection  equation.  With  a  bit  of  manipulation,  (16)  can  be 
transformed  into  a  system  of  ordinary  differential  equations  given  by: 


(18) 


In  order  to  solve  (17),  a  grid  must  first  be  formed  for  the  solution  domain.  Then,  the  state 
vector  that  corresponds  to  a  particular  grid  point  and  the  conditional  density  that  goes 
along  with  it  are  propagated  at  the  same  time.  Since  these  are  ordinary  differential 
equations,  this  can  simply  be  done  by  using  an  ordinary  differential  equation  solver,  such 
as  the  4th-order  Runge  Kutta  integrator.  [5] 
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By  combining  the  two  special  cases,  one  can  solve  the  general  FPE  used  for  this 
nonlinear  filtering  problem.  This  is  done  by  requiring  the  process  noise  spectral  density 
matrix  to  be  discrete  time.  In  essence,  the  situation  being  viewed  here  is  the  case  where 
process  noise  is  “injected”  into  the  system  at  discrete  time  intervals,  which  coincide  with 
the  measurement  times.  In  mathematical  form,  this  means: 

a=Z os) 

k 

Where  8(t-tk )  is  the  Dirac  Delta  function.  By  letting  Qt  =  ^Qk8(t -tk)  one  can  say 

k 

that  immediately  after  the  measurement  step,  the  FPE  reduces  to  the  diffusion  equation 
(15),  and  between  measurements  the  FPE  reduces  to  the  convection  equation  (16).  Since 
the  solutions  for  both  special  cases  have  been  described,  the  FPE  for  this  nonlinear 
fdtering  algorithm  can  be  solved.  [5] 


This  technique  of  using  numerical  solutions  of  the  Fokker-Planck  equation  to 
solve  the  nonlinear  tracking  problem  can  be  simply  viewed  in  the  following  flow  chart: 


http://ieeexplore.ieee.Org/iel5/7/17885/00826335.pdf 

[6] 
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One  of  the  other  areas  that  I  looked  at  this  summer  was  the  TENET  (Techniques 
for  Nonlinear  Estimation  of  Tracks)  software  developed  by  Musick  and  Greenewald.  [7] 
This  piece  of  software  implemented  two  tracking  algorithms.  One  was  a  particle  filter 
method,  while  the  other  used  an  ADI  (Alternating  Direction  Implicit)  finite  difference 
method  to  solve  the  FPE.  My  focus  was  more  on  the  ADI  finite  difference  method,  since 
it  used  a  numerical  approximation  of  the  FPE  to  solve  the  tracking  problem.  The  FPE  it 
solved,  defined  by  sub-operators,  is  given  in  the  following  form: 


dp 
d  t 


Where: 


An  alternating  direction  implicit  scheme  with  finite  difference  methods  is  used  to  solve 
the  FPE  in  this  software.  After  running  the  software  with  one  Monte  Carlo  run,  my 
results  are  given  in  the  following  diagrams: 
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FPE  Method 


PF  Method 


**Note:  Y-Axis  scales  are  different 

[7] 


FPE  Method 

**  Note:  Y-Axis  scales  are  different 
[7] 


PF  Method 


By  comparing  the  error  in  the  particle  method  against  that  of  the  FPE  based  method  (ADI 


method),  the  following  table  was  developed: 
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Error 

FPE 

PF 

Position  Error 

.476  meters 

.290  meters 

Velocity  Error 

.748  m/s 

.312  m/s 

As  one  can  easily  see,  the  particle  filter  outperformed  the  FPE  based  method  in  both 
position  and  velocity  error.  This  can  be  explained  by  the  scenario  tested  in  the  TENET 
simulator.  In  this  simulation,  a  single  dim  target  was  tracked  in  2D  space  with  a  high 
degree  of  nonlinearity.  The  a  high  degree  on  nonlinearity  and  presence  of  only  a  single 
target  creates  a  situation  favoring  the  algorithm  in  the  particle  filter  over  that  of  the  FPE 
based  method.  In  the  future,  more  work  could  be  done  with  the  TENET  software, 
especially  in  looking  into  defining  situations  where  the  FPE  based  method  outperforms 
the  particle  filtering  method,  computation  of  flops  in  order  to  compare  which  algorithm  is 
more  efficient,  and  comparing  the  performance  of  each  algorithm  against  the  resources 
used  up  in  the  running  of  the  algorithm. 

Designing  nonlinear  tracking  filters  has  been  a  very  difficult  task  for  engineers. 
The  most  direct  way  to  solve  the  problem  seems  to  be  by  using  numerical  solutions  to 
solve  the  Fokker-Planck  equation  for  the  conditional  probability  density,  but  this  quickly 
becomes  unreasonable  due  to  the  fact  that  the  solution  suffers  from  the  “curse  of 
dimensionality.”  Daum  developed  a  filter  which  uses  the  Fokker-Planck  equation  for 
probability  density  in  combination  with  an  exponential  function  and  sufficient  statistics  to 
solve  the  nonlinear  tracking  problem,  but  his  conditions  and  assumptions  are  still  fairly 


85 


limiting  to  the  real  world  application  of  the  theory.  Hopefully,  in  the  future,  a  nonlinear 
tracking  filter  will  be  developed  that  defeats  the  “curse  of  dimensionality,”  and  is  able  to 
be  applied  to  a  wider  range  of  real  world  scenarios. 
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List  of  Acronyms 


AFRL/IF 

Air  Force  Research  Laboratory  /  Information  Directorate 

ADI 

Alternating  Direction  Implicit 

AIAA 

American  Institute  of  Aeronautics  and  Astronautics 

ANTS 

Adaptive  Nonlinear  Tracking  System 

BRAT 

Baseline  Road  Assisted  Tracker 

BYU 

Brigham  Young  University 

C 

constant 

CAFE 

Complementary  Advanced  Fusion  Exploration 

C2ISR 

Command  Control  Intelligence  Surveillance  &  Reconnaissance 

CT 

Constant  Turn 

CV 

Constant  Velocity 
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DoN 

Degree  of  Nonlinearity 

EKF 

Extended  Kalman  Filter 

EKF-IMM 

Extended  Kalman  Filter-  Interacting  Multiple  Model 
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Electronic  Intelligence 
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Electronic  Support  Measures 
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Fokker  Planck  Equation 
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Fusion  Techniques  and  Nonlinear  Filtering 
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Gaussian  Mixture  Model 
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IEE 
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IHPFAT 
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Interacting  Multiple  Model  Extended  Kalman  Filter 
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Minimum  Mean  Square  Error 
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OOSP 
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